University of Birmingham Hierarchical modelling of elastic behaviour of human enamel based on synchrotron diffraction characterisation

Human enamel is a hierarchical mineralized tissue with a two-level composite structure. Few studies have focused on the structure–mechanical property relationship and its link to the multi-scale architecture of human enamel, whereby the response to mechanical loading is affected not only by the rod distribution at micro-scale, but also strongly inﬂuenced by the mineral crystallite shape, and spatial arrangement and orientation. In this study, two complementary synchrotron X-ray diffraction techniques, wide and small angle X-ray scattering (WAXS/SAXS) were used to obtain multi-scale quantitative information about the structure and deformation response of human enamel to in situ uniaxial compressive loading. The apparent modulus was determined linking the external load and the internal strain in hydroxyapatite (HAp) crystallites. An improved multi-scale Eshelby model is proposed taking into account the two-level hierarchical structure of enamel. This framework has been used to analyse the experimental data for the elastic lattice strain evolution within the HAp crystals. The achieved agreement between the model prediction and experiment along the loading direction validates the model and suggests that the new multi-scale approach reasonably captures the structure–property relationship for the human enamel. The ability of the model to predict multi-directional strain components is also evaluated by comparison with the measurements. The results are useful for understanding the intricate relationship between the hierarchical structure and the mechanical properties of enamel, and for making predictions of the effect of structural alterations that may occur due to the disease or treatment on the performance of dental tissues and their artiﬁcial replacements.


Introduction
Enamel, a highly mineralized substance, is a hard and brittle outer layer above the amelo-dentinal junction (ADJ) that covers the crown portion of the tooth, serving as an important stiff and wear-resistant part of teeth. It is a biological composite with a complex hierarchical structure, mainly composed of high content of inorganic biological hydroxyapatite (HAp) crystals (approximately 25-30 nm of thickness). Enamel consists of 80-90% inorganic material by volume (around 96% by weight) (Bechtle et al., 2012) and relatively low content of organic globular protein (Ang et al., 2012;He, 2008). At the macro-scale level, the enamel can be seen as a continuum, while at the micro-scale some notable features are present, in particular, aligned long prisms (or rods) with a specific ''keyhole'' shape, with the top oriented toward the crown of the tooth (Macho et al., 2003). On the nano-scale level, enamel is modelled as a composite with ribbon-like HAp particles organized and bundled together by globular protein (Kerebel et al., 1979).
Understanding the effect of microstructural alterations of enamel on the performance of teeth requires the knowledge of how the mechanical properties are related to the complex hierarchical structure. Over half a century, research has been carried out on the mechanical properties of enamel, and on its macroand micro-structure (He, 2008;Nakamura et al., 2011). However, few studies have focused on the influence of the nano-scale structure (Bechtle et al., 2012;Habelitz et al., 2001), where the macroscopic mechanical response can be considered to be a function of crystal shape and orientation of the mineral phase (Sui et al., 2013). Investigations at the nano-scale require the use of advanced techniques and systematic models to establish a firm basis for understanding the hierarchical structure-property relationships.
One suitable technique for the study of nano-structure is small angle X-ray scattering (SAXS), an advanced non-destructive technique widely used to reveal quantitative information about the orientation and degree of alignment in crystalline and amorphous materials (Fratzl et al., 1996). Another X-ray technique, wide angle X-ray scattering (WAXS), is widely used to quantify the internal strain in crystal lattices, including in response to external loading, e.g. in situ compression or tension (Young et al., 2007a,b). SAXS/ WAXS have been applied only recently to the study of mineralized biological composites, such as bones and bovine teeth Stock, 2005, 2010;Deymier-Black et al., 2010;Singhal et al., 2012). Deymier-Black et al. (Deymier-Black et al., 2010) determined the longitudinal apparent modulus of hydroxyapatite (HAp) in bovine dentine using synchrotron based WAXS, while strain distribution across the ADJ in bovine teeth was investigated by Almer and Stock . None of these studies were performed on human enamel, in which different particularities of its morphology are expected to result in differences in the mechanical properties (Nogueira et al., 2011). Moreover, very few studies devoted to human enamel have been published (Fujisaki et al., 2012), and early studies did not take into account the nanoparticle shape, size and orientation distribution (Bechtle et al., 2012). Therefore, a deep understanding of the relationship between the nano-scale structure and the macroscopic mechanical behaviour of human enamel is still lacking.
Many models of composite mechanical behaviour have been proposed, notably the Voigt and Reuss bounds, as well as the Jones (Jones, 1999) and BW models (staggered microstructural model) ( Bar-On and Wagner, 2012). Besides these, an improved multi-level model, based on the original Eshelby inclusion model (Takao and Taya, 1987;Withers et al., 1989) has been established and successfully applied in the evaluation of the mechanical response of human dentine by capturing the relationship between the nanoscale structure and macroscopic loading (Sui et al., 2013).
In this study, the in situ synchrotron technique, simultaneous SAXS/WAXS, was used to measure the elastic crystallite strain (WAXS), the alterations in crystal orientation and the degree of alignment of nano-scale HAp crystals (SAXS) in human enamel subjected to an externally applied uniaxial compressive loading along the approximately longitudinal direction with respect to the preferential rod orientation. The multi-scale Eshelby model was modified and extended to capture the hierarchical structure and properties of the human enamel. The model predictions were verified by comparison with experimental data. The ability of the model to capture the observed response is discussed.

Sample preparation
A freshly extracted normal human third molar (ethical approval obtained from the National Research Ethics Committee; NHS-REC reference 09.H0405.33/Consortium R&D No. 1465) was washed and cleaned in distilled water to eliminate residues so that the possibility of contamination or other chemical effect was excluded. The samples were kept for a maximum of 7 days in distilled water in a commercial fridge at 4°C until the experiment was performed, implying relatively low diffusion rates that impede or exclude any possible superficial demineralization effect of water. It is also important to highlight that the technique we employed (X-ray diffraction in transmission) is a bulk, not a surface technique. As a consequence, even if small superficial changes were to arise, their contribution to the overall scattered signal would be negligible. Thus, the composition of the inorganic phase in the dental tissue was well-represented by pure organically-derived HAp. The crown was cut off below the ADJ using a low speed diamond saw (Isomet Buehler Ltd., Lake Bluff, Illinois, USA) and further cut into smaller cubes, and a series of polishing papers were used to produce the final 2 Â 2 Â 2 mm cube of enamel. Regarding the rod direction of the enamel sample, the cut in the prepared sample was close to the region perpendicular to the occlusion surface. Therefore, the predominant rod direction was assumed to be in the longitudinal direction of the sample.

Micro-CT protocol and data processing
To determine the measuring positions and loading cross-section for the SAXS/WAXS experiments, a micro-CT scan was carried out at 1.9 lm resolution using 40 kV voltage, 120 lA current and a 0.5 mm Aluminium filter (SkyScan 1172, Bruker microCT, Kontich, Belgium). The resulting slices were reconstructed using the Sky-Scan NRECON software package, and subsequent 3-D planning models were created using Fiji imaging software (Eliceiri et al., 2012).

In situ scattering measurements
The experiment was performed on the B16 beamline at Diamond Light Source (DLS, Oxford, UK). A schematic diagram of the experimental set-up is shown in Fig. 1a. The sample of human enamel was slowly deformed in compression along the x-direction in the laboratory coordinates ( Fig. 1a) at the displacement rate of 0.2 mm/min until failure, using a remotely operated and monitored compression rig (Deben, Suffolk, UK), with a 5kN calibrated load cell. The rig was equipped with custom-made jaws, allowing a high-energy transmission X-ray setup to be used.
The monochromatic X-ray beam at the photon energy of 17.99 keV was used and collimated to the spot size of 0.5 Â 0.5 mm. The beam was incident at the sample in the direction perpendicular to its loading direction. Two separate WAXS and SAXS detectors were alternately setup to collect the patterns at consecutive loading increments downstream of the beam. The global Z-axis is along the beam direction; global X-axis corresponds to the horizontal (loading) direction, Y to the vertical. The X-Y plane is parallel to the detector plane. WAXS diffraction patterns were recorded using a Photonic Science Image Star 9000 detector (Photonic Science Ltd., UK) placed at a sample-to-camera distance of 128.72 mm. Further downstream, a Pilatus 300 K detector (Dectris, Baden, Switzerland) was positioned at a distance of 4358.47 mm to collect the SAXS patterns. In order to record both the WAXS and SAXS patterns at each scanning location, the WAXS detector was translated laterally to expose the SAXS detector after each WAXS exposure. A lightly compacted disk of standard silicon powder and a dry chicken collagen sample inserted close to the sample position were used as calibration standards and to determine the sample-to-detector distance with the required precision (Calibration, 2013).

Scattering data analysis
2.4.1. WAXS data analysis WAXS data can be interpreted in terms of the shift of the diffraction peak obtained from a cluster of HAp crystals, so that the average micro-strain (lattice strain) in the crystals can be deduced (Young et al., 2007a,b). The typical WAXS pattern of HAp is shown in Fig. 2a (only the (0 0 2) peak is selected for interpretation). At each compressive load step, the beam was scanned between the loading platens, and only the results from the middle were selected for WAXS interpretation to guarantee that the same location was studied during compressive deformation. The elastic lattice strain of the HAp phase was computed by calculating the changes in the inter-planar spacing between the crystal planes (Korsunsky et al., 2011): where d 002 is the deformed d-spacing and d 0 002 is the reference strain-free value. In detail, 2D diffraction images were firstly preprocessed using Fit2D (Hammersley, 1997). The (0 0 2) peak of interest from each pattern was ''caked'' (a professional jargon 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 1D radial pattern) within the range of 20°around the loading direction (Fig. 2a). The strain measurement needs high quality of peaks in order to determine the shift of the peak centre by Gaussian fitting. The smaller the angular width of the bin is, the lower (or noisier) the integrated (averaged) intensity will be. Since the crystalline perfection in biomaterials is not as good as that in the pure crystal, the binning angle of 20°is often used in order to guarantee the pronounced peak to fit. Furthermore, each pattern was ''caked'' with the same width of 20°in the direction of À15°, 0°, 15°, 30°and 45°of the (0 0 2) peak (see Fig. 2a). The normal strain component along the centre direction of each ''cake slice'' represents the strain component in the corresponding orientation. Subsequently the 1D radial plot of each individual (0 0 2) peak within each sector was fitted with a Gaussian function to obtain the peak centre position. The sample under strain-free condition (without any load) was characterised by WAXS and used as the reference state for strain measurement. As the load increased, the shift of the peak centre position with respect to the strain-free reference point allowed the calculation of the HAp elastic lattice strain. In addition, the structural orientation angle was determined from the strain-free sample by azimuthal-radial ''caking'' of the (0 0 2) peak over the range of available peak (À30°$ 60°), and fitting the azimuthal centre position of the pronounced peaks of I WAXS (u)$u_ WAXS plot ( Fig. 2c) (Al-Jawad et al., 2007).

SAXS data analysis
For SAXS data analysis, the pattern from the strain-free enamel sample was taken as a reference. Due to the dense packing of crystals in enamel, it is the electron density change occurring in the gaps between crystalline particles that gives rise to the scattering signal (Tanaka et al., 2010). It is also understood that the orientation of the gaps between rods roughly coincides with the orientation of the crystals within the rod (Tanaka et al., 2010). Thus, the information from gap scattering can be used to deduce the orientation and degree of alignment (percentage of aligned particles) of HAp crystals. To quantify it, the 2D SAXS patterns were processed by integrating over the entire relevant range of scattering vector q, resulting in a function I SAXS (u) of the azimuthal angle u_ SAXS ( Fig. 2d) (Tesch et al., 2001(Tesch et al., , 2003. The predominant orientation / 0 SAXS of the mineral crystals is determined by the position of the two peaks in the plot of I SAXS (u) (e.g. / 0 SAXS in Fig. 2d), which is corresponding to the short axis of SAXS pattern. Further, the degree of alignment q with respect to the predominant orientation of HAp can be calculated by the ratio of the two areas under the curve of I SAXS (u): where A unoriented is the area of constant background level accounting for the scattering from unoriented particles and A oriented depicts the total area below the curve of I SAXS (u)$u_ SAXS subtracting A unoriented . Thus, the value of q ranges from 0 to 1, with q = 0 indicating no predominant orientation within the plane of the section, while q = 1 indicates a perfect alignment of all the crystals (Tesch et al., 2001(Tesch et al., , 2003.

Geometrical assumptions
Human enamel has a hierarchical two-level composite structure, where the first level is represented by a rod with a keyholelike section and the second level by the organized and bundled HAp crystallites within each rod (Nanci and Ten Cate, 2003). Fig. 3a shows the keyhole-like microstructure of enamel, modified from Habelitz (Habelitz et al., 2001) and the distribution of the HAp crystals within the rod in 3D. Fig. 3b-d provides schematic illustration of the geometric model derived from the enamel structure, where the first-level regards the whole enamel sample as composed of aligned rods bundled by the joining medium of globular protein phase (Fig. 3b), and the second level considers the rod as a composite in detail, consisting of partially aligned HAp crystals and globular proteins.
In theoretical model, especially in describing complex systems, some reasonable simplifications need to be made in order to capture the basic mechanism. Thus the direction of rods in the first level is assumed to be uniform to minimize the number of variables. Based on the objective of this study, which is the nano-scale particle distribution effect on macro-scale elastic behaviour, WAXS and SAXS could capture the orientation distribution (i.e. texture) of HAp crystals in the second level, which in turn can be used to reflect the orientation inside actual rod in the gauge volume since the rod is composed by HAp crystallites. The partially aligned distribution of HAp crystals is based on the observation using scanning transmission electron microscopy (STEM), where a high percentage of HAp crystals have a perfect alignment while the other crystals are more randomly distributed (Huang et al., 2010). Both levels are non-dilute systems consisting of a number of inhomogeneous inclusions. For simplicity, both rod and HAp crystals are assumed to be of ribbon-like shape (Fig. 3b-d).
In a previous study of dentine, the penny-shape inclusion was used to simulate the platelet-like structure of HAp crystals in human dentine (Sui et al., 2013). As for the ribbon-like shaped HAp crystals in human enamel, a ribbon-shaped inclusion is used to simulate an individual HAp crystal. In the Eshelby approach, the cylinder is approximated by a prolate spheroid described by the three dimensions, a 1 , a 2 and a 3 in Fig. 3e. Normally a 1 = a 2 ( a 3 , i.e. the cross section is a circle (Mura, 1987). The crystal c-axis (corresponding to the (0 0 2) peak) of ribbon-like shaped HAp crystals is normal to the cross section of the cylinder, i.e. along a 3 direction (Huang et al., 2010). In the next section, the model is briefly introduced, while a detailed derivation is shown in Appendix.

First-level model: multiple aligned rod inclusions within 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 (Clyne and Withers, 1995), the stress in the rod can be expressed as In the above expression, f 1 is the volume fraction of rods in the enamel, S 1 is the Eshelby tensor for a cylinder that approximates the rod shape, C matrix , C rod are the stiffness tensors of the globular protein and rod respectively, and T is the tensor transformation (rotation) matrix that depends on the Euler angles (see Appendix) giving the orientation of an individual rod with respect to the laboratory coordinate system. In the present model, the laboratory coordinates were fixed, and it was assumed that the rods were all aligned along the loading direction, thus T was constant. The rod stiffness remains to be determined from the second level model.

Second-level model: partially aligned HAp inclusions within each individual rod
The purpose of the second-level model of enamel, which describes partially aligned HAp inclusions within one rod, is to establish the relationship between the rod inclusion stress and the average strain in the HAp crystals in rod hei HAp (r rod = K 1 hei HAp ). The measured crystal strain corresponds to the mean strain value for all the crystals within the considered gauge volume (Chou et al., 2012). Due to the partial alignment of HAp crystals within the rod, the real apparent stiffness K 1 is to be given by the values bounded by the two extreme cases, namely that of fully random distribution and that of perfect alignment.

Multiple perfectly aligned HAp crystals
Supposing all the crystals are perfectly aligned and the alignment direction is described by the transformation matrix T, the relationship between the average HAp strain and the rod inclusion stress obtained from the first-level model can be established as follows: where hCi HAp and hS 2 i are the average stiffness and Eshelby tensor. Since all the crystals are perfectly aligned, hCi HAp and hS 2 i can be represented by the single crystal values hCi HAp = C HAp , hS 2 i = S 2 . Note that a single HAp crystal in the second level can be regarded as a single inhomogeneity with ribbon-like shape, thus the Eshelby tensor for a single HAp should be the same as the rod inclusion (S 2 = S 1 ) and the stiffness of a single HAp has been reported to be transversely isotropic (Ochsner and Ahmed, 2011). In Eq. (4), different orientation angles (different T) would lead to different values of hei HAp aligned .

Multiple randomly distributed HAp crystals
Supposing all the crystals are randomly distributed the relationship between average HAp crystal strains and rod stress is independent of the transformation matrix in Eq. (4). Thus the relationship becomes matrix r rod or expressed more simply r rod ¼ K random hei HAp random Different from perfectly aligned crystals, the average value of hS 2 i can no longer be the value of a single crystal, but should be obtained by the volume average value of all the crystals. However, as an alternative, the averaging effect can be captured by using the single crystal relationship as Eq. (4) and averaging over all the values with different Euler angles (see Appendix, averaging the results obtained from each single crystal relationship as Eq. (4) over all possible orientations).

The combination of random and aligned HAp inclusions
In order to model the observed preferred orientation of HAp inclusions (Fig. 3d), we represent it as the combination of a random distribution with a volume fraction f aligned of perfectly aligned particles. In this case the overall apparent stiffness K HAp partial aligned is given by the rule of mixture between K random and K aligned : where (1 À f aligned ) now represents the volume fraction of the randomly distributed, and f aligned corresponds to the degree of alignment of crystals revealed by SAXS. Thus, the relationship between rod inclusion stress and average internal HAp lattice strain becomes

Determination of the rod stiffness
In Eq. (3), the stiffness of the rod C rod still remains to be determined. Since the rod is regarded as a composite consisting of globular protein and mineral HAp crystals, C rod can also be determined using the Eshelby inclusion model, taking into account the volume fraction (f 2 ), average Eshelby tensor (hS 2 i) and average stiffness (hC HAp i) of HAp crystals in the second level. To simplify the determination, it is assumed that all HAp crystals in the rods are aligned with an angle with respect to the longitudinal direction of rod (thus hCi HAp = C HAp , hS 2 i = S 2 ). Such orientation is the preferred orientation revealed by SAXS/WAXS and is reflected in the model by a unique transformation matrix T HAp . The expression of the stiffness of rod is given here without detailed derivation whereS 2 ¼ T T HAp S 2 T ÀT 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. Further, the overall relationship between the average HAp strain and the externally applied stress can be established by combining Eqs. (3) and (7).
4. Experiment result and model evaluation

Nano-scale HAp distribution and mechanical response of enamel
The loading areas of the samples were accurately determined by micro-CT measurement, which is 3.55 mm 2 . Fig. 2a shows a WAXS pattern of enamel consisting of a system of Debye-Scherrer rings (peaks). Since the enamel is textured, only limited range of rings can be captured. The apparent radial shifts of the (0 0 2) peak in the WAXS pattern were measured under uniaxial compressive loading applied in longitudinal direction with respect to the rod direction. The preferential orientation of HAp crystals obtained by WAXS pattern is shown to be roughly perpendicular to the arc of (0 0 2) peak and the detailed value is determined from the stress-free plot of I WAXS (u)$u (Fig. 2c) with the Gaussian fit (red line).
The SAXS pattern as shown in Fig. 2b contains the information of orientation and degree of alignment of the gaps between the HAp crystals. As mentioned above, the information about HAp crystallites orientation distribution can be deduced by gap scattering. The preferential orientation shown in Fig. 2b is roughly along the short axis of the elliptical pattern and the detailed value is also determined from the stress-free plot of I SAXS (u)/u ( Fig. 2d) with the Gaussian fit (red line). The detailed values of the orientation and degree of alignment obtained by WAXS and SAXS patterns are listed in Table 1, from which these parameters were used in the Eshelby model evaluation.
The shifts of (0 0 2) peak from WAXS along x-direction were used to obtain the elastic lattice strain values. Fig. 4 shows the experimental results of the applied stress vs. HAp lattice strain of sample along the loading direction, indicating a linear increase as expected. Two results are given respectively using the preferred orientation angles obtained by SAXS and WAXS but keeping the other parameters the same. The ratio of the uniaxial stress and the average HAp lattice strain gives the apparent modulus (Powers and Farah, 1975), which is listed in Table 1. Meanwhile, the residual (initial) strain was found to be quite small and can therefore be neglected. The shifts of (0 0 2) peak along other directions were also measured, e.g. by caking each pattern with the width of 20°along the directions of À15°, 0°, 15°, 30°and 45°to determine the azimuthal variation of the normal strain component. The result is shown as an azimuthal plot in Fig. 5a and b, where 0°represents the loading direction. Note that the normal strain component undergoes a transition from negative and to positive at around 25°.

Evaluation and testing of multi-scale Eshelby model
In the model, the material properties and other parameters used were taken from the literature, and refined by fitting with the experimental data. The average mineral concentration (HAp volume fraction) has been reported to be $95% at each level of enamel (Bechtle et al., 2012). In detail, at the first level, 95% represents the volume fraction of rods within the globular proteins, while at the second level, 95% represents the volume fraction of mineral within each rod. This means that the overall volume fraction of mineral in the enamel is 95% Â 95% $90%. In general, Young's modulus of 1 GPa for globular protein is found in the literature, without taking into account the viscoelasticity or viscoplasticity (Huo, 2005;Qin and Swain, 2004). Polycrystalline HAp is considered to be transversely isotropic with five independent elastic constants (Ochsner and Ahmed, 2011). To represent the shape of the rod and of the HAp crystallites for each level, the Eshelby tensor for the cylinder (prolate spheroid) was used. The elliptical semi-axes a 1 and a 2 within the transverse cross-section were assumed to be the same, (a 1 /a 2 = 1), but different from a 3 . The apparent modulus K was calculated based on the different preferred orientation angles of the HAp crystals obtained by both SAXS and WAXS. All the parameters refined to obtain the best fit are listed in Table 1, as well as the reference values from literature. A comparison of the stress/strain curve along the loading direction between the model prediction and the experiment is plotted in Fig. 4, where the Voigt and Reuss bounds are also indicated. As expected, it is found that the modified Eshelby model prediction lies between the two bounds and gives a satisfactory agreement with the experimental data. Meanwhile, the comparison of normal strain variation with the azimuthal angle is presented in Fig. 5a and b. Fig. 5a gives the model prediction results using different preferred orientation angles of HAp crystals obtained by both SAXS and WAXS. It is apparent that the model prediction of the transverse tensile strain that arises under compression (the Poisson effect) falls short of the observed strain. A satisfactory agreement can be achieved by adding a tensile transverse strain component (along the global Y-axis in Fig. 1) perpendicular to the rod and loading direction (the global X-axis in Fig. 1). The updated result is illustrated in Fig. 5b. This remarkable effect deserves further detailed analysis and discussion.

Discussion
The Eshelby model can be constructed to simulate the hierarchical structure. It has the strong physical significance in that it captures the interaction between two different phases (inhomogeneity and matrix), compared with other simple models, like Voigt and Reuss. For the interface, the Eshelby model assumes fully coherent interface, in line with the prevailing view in the literature. The main shortcoming of the shear-lag model (Ji and Gao, 2004) is that it does not consider the influence of the nanoparticle distribution within the enamel. This is the key aspect that Eshelby model allows deeper understanding of the relationship between the nano-scale structure and mechanical behaviour, even though the uniform Eshelby tensor is used to describe the monodisperse distribution of crystal size for the simplifications when developing a theoretical model. Furthermore, the key aspect of interest in this study lies in the apparent modulus of enamel, i.e. the relationship between the strain in the nano-crystallites and the macroscopic material loading. The composite effect cannot be ignored. The woven structure (Palmer et al., 2008) of blocking effect is expected to be effective only during tensile loading rather than the uniaxial  compressive loading. It is considered particularly when investigating the fracture mechanism of enamel rather than the elastic deformation.

Refined parameters by Eshelby model
The refined parameters of the multi-scale Eshelby model are listed in Table 1. Most of the values reported in the literature were found to give satisfactory results. The key parameters varied in the analysis were the volume fraction (the volume fraction of rods in the enamel and of HAp crystallites within rods were assumed to be the same) and the elastic constants of the HAp crystal. Among these parameters, the elastic constants have the most significant influence on the result. In the optimization process, the volume fraction was first fixed at the reported value (95%), and the elastic constants, especially the Young's modulus, were refined. Subsequently, other parameters (e.g. the globular protein modulus) were refined as well, although it was found that they have a minor effect on the apparent result. This leads to the conclusion that the elastic modulus of soft globular protein phase, whether it be taken as isotropic or anisotropic, does not exert a large effect on the validation within the scope of model-experiment matching considered here. In addition, the refinement process also helps in the identification of nano-scale parameters, which may be hard to determine directly from the experiments. In addition, the model can be applied to diseased tissues as well, in which case, the model needs to be adjusted correspondingly, given different initial parameters to be refined.

Residual strain
During the preparation process using a low speed diamond saw and polishing papers, as even during the natural growth, a thin layer of residual strain may be induced at the sample surface. However, as shown in Fig. 4, the pre-existing residual strain is very small. Since only the linear elastic response of enamel is considered in the experiment (reflected in the linearity of the experimental stress-strain curve), the presence of initial residual strain only amounts to an offset that does not affect the parameters like the apparent modulus. Hence, the low level residual strain was neglected in the present analysis.

Gap scattering of SAXS
The good agreement between the apparent modulus K results calculated using different HAp preferential orientations determined by SAXS and WAXS provides strong evidence and validates the argument that the SAXS pattern arising from gap scattering (Tanaka et al., 2010) can be used to deduce the HAp crystallite orientation distribution, i.e. the gaps are almost parallel to HAp crystals inside the rod. Therefore, we report direct confirmation that SAXS data provides beneficial complimentary information for determining the crystal orientation distribution in enamel and other mineralized tissues.

Normal strain components variation
The normal strain components variation of the HAp crystals with respect to different azimuthal angles (0-90°) is shown in  Fig. 5b between the experimental data and the modified model incorporating transverse strain) indicates that the deformation state in the sample is likely to be multi-axial, with an additional tensile strain component existing along the transverse direction. An explanation for this observation must lie in the fact that a disproportionately large transverse tensile strain is developed within the sample due to the interaction between the material structure and the loading arrangement. One possible mechanism to explain this is the barrelling effect (mid-section expansion) in the enamel sample: the external uniaxial compression causes local locking of the enamel to the platen surface accompanied by the expansion of the sample perpendicular to the loading direction in the midsection (where the measurements were performed). In addition, micro-cracks smaller than 2 lm (below the resolution of the micro-CT scan) pre-existing in the enamel may result in the debonding between rods (Bajaj and Arola, 2009), and thus modify the transverse strain. Such mechanisms are not captured by the multi-scale Eshelby model and need to be introduced externally. The validation of the precise mechanism needs direct observation using advanced ultra-high resolution (sub-micron) imaging and analysis techniques, e.g. digital volume correlation.

HAp crystals orientation distribution effects
The combination of random and perfectly aligned HAp inclusions (that can be referred to as ''partially aligned'') was described in the previous section (Fig. 5). We note that this is an approximation, since in practice the orientation distribution of HAp crystals is a continuous function of the angle between the c-axis of HAp crystals and a chosen direction (the global X-axis) in the laboratory coordinate system. The WAXS intensity variation with the azimuthal angle plotted in Fig. 2c allows us to extract the volume fraction of HAp crystals with different orientations. Based on this, continuous orientation distribution model can be carried out based on the probability density function of finding crystals of given orientation. The result obtained using this approach is shown in Fig. 6. The comparison between this result and that obtained earlier using the combination of random and perfectly aligned crystallites indicates that the difference between the two approaches is negligible in terms of the prediction of the apparent modulus. Thus, the angle of the preferred orientation has a dominant effect on the mechanical response.

The effects of preferred orientation of HAp crystals
The effect of the nano-scale structure (the HAp crystallite distribution) on the macroscopic mechanical response was further investigated by changing the preferential crystal orientations (changing the transformation matrix in Eq. (5)). A 3D model of perfectly aligned crystals inside a rod is established (Fig. 7a) with the angle u describing the rotation of the alignment direction around the global Z-axis. The local coordinate system with respect to the apatite has been designated x,y,z in Fig. 7a. When all HAp crystals are aligned along the global X-axis, u equals to 0°. By changing the relative alignment direction, the variation of K aligned with respect to the loading direction can be calculated (Fig. 7b). From Fig. 7b, the corresponding results using the real orientation angles found in the experiment (174°from SAXS and 14°from WAXS) are found to be Kaligned SAXS ¼ 138:1 GPa and Kaligned WAXS ¼ 131:8 GPa, i.e. closely similar values. Meanwhile, due to the high degree of alignment of HAp crystals in the enamel, the value of the overall apparent modulus K HAp partialaligned lies closer to K aligned rather than K random . The enamel displays strong microscopic elastic anisotropy. It is interesting to note that the orientation with the largest stiffness is found, as expected, around 0°with respect to the loading direction. However, the most compliant orientation observed is not at 90°(perpendicular to the loading direction), but rather around 50°or 130°. This is due to the transversely isotropic stiffness of HAp crystals.

Conclusions
In this study, the longitudinal apparent modulus of human enamel was measured during in situ elastic compression by the combination of synchrotron WAXS and SAXS. This is the first combined SAXS/WAXS study on the nano-scale structure and its influence on the macroscopic mechanical behaviour of human enamel. The coherent results obtained for both WAXS and SAXS indicates that the SAXS pattern arising from the inter-mineral gaps can be used to reflect the HAp crystal orientation. This provided access to the information on both the structural and mechanical aspects of the sample and allowed us to make further progress compared with the previous studies which only used WAXS (Almer and Stock, 2005;Fujisaki et al., 2012). As an improvement of the earlier proposed composite model ( Bar-On and Wagner, 2012;Jones, 1999), a multi-scale Eshelby inclusion model was established to estimate the elastic properties of enamel, considering it as a two-level The beam direction is along the global Z-axis and the alignment here represents the angle between the local x-axis and global X-axis (initially the aligned angle u = 0°, i.e. the local x-axis of a ribbon-like shape is initially along the global X-axis) and (b) by changing the alignment angle (0-180°), the average strain of the crystals along the loading direction can be obtained by the multi-scale Eshelby model, which then indicates that K aligned varies with respect to the preferential alignment angle.
composite. Good agreement between the measured and calculated lattice strains along the loading direction validates the model. With respect to the normal strain component in a general azimuthal direction, the underestimation of the tensile strain by the model required the introduction of additional transverse internal strain. The modified model result demonstrates that the deformation state in the enamel sample may not have been fully uniaxial and may have been caused by the barrelling effect or crystal debonding.
This systematic experimental and modelling approach reported here is able to capture the complete picture of the multi-scale structure, material elastic properties and its evolution under loading. The parameter refinement and validation approach adopted in the present study offers an important alternative route to the identification of nano-scale parameters. Other established experimental characterisation techniques, such as nanoindentation and microscopy, are confined to the sample surface, making overall bulk parameter identification difficult. Combining the present results on enamel with previously published data on human dentine, an improved and comprehensive understanding of the multi-scale structural-mechanical properties within human dental tissues can be given. Besides the implications for the characterisation of the structure-property relationship of other hierarchical biomaterials, this knowledge is essential for developing better prosthetic materials and dental fillings, and could also shed light on the mechanical property evolution associated with the multi-scale structural changes within human teeth due to disease and treatment.
e HAp where T is the orientation matrix described by three Euler angles (h, u, w). The partial alignment is understood to be a combination of perfect alignment and random distribution. The expression for the average strain in multiple perfectly aligned HAp crystals is similar to Eq. (A5), as shown in Eq. (4), while that for the average strain in multiple randomly distributed HAp crystals needs to be obtained by the volume average method, which is introduced here. The purpose of using volume average method is to avoid the complex calculation of average Eshelby tensor and average stiffness of randomly distributed HAp crystals. In a random distribution, each crystal follows the relationship of Eq. (A5) by an individual transformation matrix T. Crystals can have any possible orientation in the space, thus the volume average method is to calculate the mean strain value over all the crystals in the rod by averaging over all possible orientations based on Eq. (A5). e HAp