3D Maps of Mineral Composition and Hydroxyapatite Orientation in Fossil Bone Samples Obtained by X-ray Diffraction Computed Tomography

Whether hydroxyapatite (HA) orientation in fossilised bone samples can be non-destructively retrieved and used to determine the arrangement of the bone matrix and the location of muscle attachments (entheses), is a question of high relevance to palaeontology, as it facilitates a detailed understanding of the (micro-)anatomy of extinct species with no damage to the precious fossil specimens. Here, we report studies of two fossil bone samples, specifically the tibia of a 300-million-year-old tetrapod, Discosauriscus austriacus, and the humerus of a 370-million-year-old lobe-finned fish, Eusthenopteron foordi, using XRD-CT – a combination of X-ray diffraction (XRD) and computed tomography (CT). Reconstructed 3D images showing the spatial mineral distributions and the local orientation of HA were obtained. For Discosauriscus austriacus, details of the muscle attachments could be discerned. For Eusthenopteron foordi, the gross details of the preferred orientation of HA were deduced using three tomographic datasets obtained with orthogonally oriented rotation axes. For both samples, the HA in the bone matrix exhibited preferred orientation, with the unit cell c-axis of the HA crystallites tending to be parallel with the bone surface. In summary, we have demonstrated that XRD-CT combined with an intuitive reconstruction procedure is becoming a powerful tool for studying palaeontological samples.


S1. Experimental details
The XRD-CT measurements were carried out at European Synchrotron Radiation Facility (ESRF) at beamline ID15A, using a monochromatic beam of energy 86.6 keV, corresponding to a wavelength of λ = 0.143 Å. Compound refractive lenses and slits were used to focus the beam to a square spot (see below for dimensions) on the sample. A diode in the beam path upstream of the sample was used to correct for variations in the incoming photon flux. The sample was mounted on a stage allowing translations in the x and y directions along with rotation of the sample by ω about the tomographic y-axis. A beam stop was placed in front of the detector to prevent the intense direct beam from damaging the detector. The area detector was a Perkin-Elmer XRD 1621 N ES Series, operated in 2 × 2 binning mode to achieve a frame rate of 20 Hz, with 1024 × 1024 pixels and an effective pixel size of 400 m × 400 m. The time-varying dark current of the detector was accounted for by recording dark images. Slight misalignments in the setup were numerically corrected using the scattering signal from a CeO2 powder reference. Corrections for polarization and Lorentz factor were found to be negligible at the small 2θ angles studied. Absorption corrections were not made, due to the qualitative nature of the analysis.
The sample of Discosauriscus austriacus (Fig. S1a) was measured with the long bone axis parallel to the tomography axis y. The beam was focused to a 20 m × 20 m spot. The sample was scanned in x in steps of Δx = 20 m and rotated about the y-axis by the angle ω in angular steps of 0.90º. A total number of 250 different x and 200 different ω were used, where ω  [0, 180º]. Only one value of y was measured, giving a cross section of the diaphysis (Fig. S1a). These parameters provided pixels of approximate size 20 m × 20 m in the reconstructed single-slice tomogram.
For the sample of Eusthenopteron foordi, three independent full sets of XRD-CT measurements were performed with approximately orthogonal orientations of the sample on the sample stage (Fig. S1b), allowing information about the preferred orientation to be extracted. For each of the three orthogonal 3D tomographic measurements, the same measurement procedure was used. The resolution was reduced by refocusing the beam to a 150 m × 150 m spot. The position in y ("sample height") was scanned across the sample in 100 steps, separated in y by Δy = 150 m. For each position y, the sample was scanned in x in steps of Δx = 150 m and rotated about the y-axis by the angle ω in angular steps of 2.22º. A total number of 98 different x and 81 different ω were used, where ω  [0, 180º]. This provided voxels of size 150 m × 150 m × 150 m in the reconstructed 3D tomograms. The 3D image registration problem of aligning the three datasets was solved by fitting three translation and three Euler angle parameters. The 2D XRD-CT measurements of the tibia of Discosauriscus austriacus resulted in acquisition of approximately 100 GB of raw data, while the 3D XRD-CT measurements of the humerus of Eusthenopteron foordi resulted in approximately 3.6 TB of raw data. To reduce the amount of data, each diffraction pattern was divided into 64 sectors, and the intensity azimuthally averaged within each sector, thereby reducing the amount of data to 0.80 TB in total for both samples, at the cost of reducing the resolution in  to Δ = 360º/64  5.6º. For efficient data processing, the intensity maximum of each Bragg peak was used to estimate the intensity contributions to the sinograms.

S2. A detailed description of XRD-CT
In an XRD-CT experiment, the observed intensity of the diffracted beam, rather than the reduction in intensity (attenuation) of the direct beam, is used to reconstruct the sample volume (or cross-sectional area). Whereas one usually uses a wide beam flooding the whole sample to record full projection images for attenuation CT, in the case of XRD-CT the sample is scanned using a narrow beam, and a 2D diffraction pattern is recorded for each position and orientation of the sample. A schematic setup of the XRD-CT experiment is shown in Fig. 1b, where farfield diffraction is recorded by a 2D-detector. The sample can be translated along the x and y axes as well as rotated about the y-axis with the sample rotation angle ω and illuminated by a monochromatic X-ray beam to measure diffraction from different parts of the sample for different directions.
The observed scattered X-ray intensity I in the far-field is proportional to the absolute square of the Fourier transform of the electron density in the sample 1 denotes the Fourier transform of the electron density. The scattering vector is defined as Q = ks -ki, where ki denotes the incoming and ks the outgoing wave vectors. The magnitude of the scattering vector is given by Q = |Q| = 4π sin(θ)/λ, where θ is half the total scattering angle 2θ and λ the photon wavelength. Because of the high energy of the radiation only scattering to small angles 2θ is observed. The reconstruction challenge is to assign -via tomographic methods -densities to the different voxels in the numerically reconstructed sample. Because diffraction patterns tend to exhibit circular features for polycrystalline samples, it is expedient to use polar coordinates, i.e., the measured intensity is a function of Imeasured = Imeasured(ω,Q,,x,y) where  is the azimuthal scattering angle, as illustrated in Fig. 1. Intensity distributions were extracted and tomograms obtained as described in detail (vide infra), using the MATLAB® built-in inverse radon function based on the filtered back projection (FBP) with a Ram-Lak filter.

S3. Data analysis
The approach explored in this study is to project each recorded diffraction pattern onto one single scalar value ("intensity descriptor") Idescriptor(,x,y), analogous to how one single transmitted intensity parameter It(,x,y) is measured in conventional CT. A requirement is that the signal should be (at least approximately) invariant under rotations, as further discussed below. Then, we hypothesize that the computationally efficient and well-established filtered back-projection (FBP) reconstruction procedure can be inherited from attenuation CT. Clearly, there are many ways of reducing a diffraction pattern to a single scalar, as we shall elaborate in the following.
"Isotropic" descriptors A particularly simple case is that of many isotropically oriented crystallites (a powder): The scattering can then be considered 1D ("powder-like"), and if choosing a particular Q = Qhkl corresponding to a given Bragg peak, this diffracted intensity integrated over  is a natural choice for the single scalar parameter. Mathematically, this can be expressed by For example, the 002HA Bragg peak of an isotropically scattering region of HA can be integrated over all ϕ for each diffraction pattern, giving a reduced dataset I 002 isotropic (ω,x,y), which can then be used to reconstruct in 3D where this crystalline phase is located inside the sample. This simplistic approach to diffraction contrast imaging is surprisingly powerful, and has even been used successfully in cases for a similar technique to XRD-CT, namely small-angle X-ray scattering computed tomography (SAXS-CT) where the scattering substance is not fully isotropic 2 . Note that for improved statistics, it is also in many cases possible to choose the summed intensity of several Bragg peaks, appropriately scaled for the unit cell structure factor, as the scalar parameter, a technique utilized in this work for reconstructing the tomograms of different minerals in the Eusthenopteron foordi sample (Fig. 4). The total scattered signal Itotal is obtained by integrating Imeasured over the whole detector, i.e. In the more interesting case of preferentially oriented crystallites, the directional (vector) properties of Q must be accounted for. With an ensemble of crystallites exhibiting a shared preferred orientation, the diffraction pattern generally changes as the sample rotates. As mentioned, this is a complicating factor for tomography, while it also opens for vectorial imaging. If the sample contains only a few large single-crystalline grains, the algorithms described here are inadequate, as the diffraction signal will change too abruptly in ω and ϕ. We note that other groups have developed complementary methods for dealing with this special case 3 . However, if the crystallites are numerous, have a sufficiently low degree of preferred orientation, and small orientation gradients, it is possible to construct intensity descriptors that enable vectorial tomographic reconstructions based on filtered back-projection. Figure S2 a) Illustration of the Ewald sphere for a polycrystalline sample with preferred orientation. The sketch illustrates the case of crystallites with a broad and continuous orientation distribution, such as the HA crystallites studied here. b) The 002HA peak exhibits pronounced preferred orientation, with a peak width Δ ~ 50º. c) Quartz Bragg peak. d) Calcite Bragg peak. The signals are qualitatively different: While the quartz 100qu Bragg peak is essentially isotropic, the calcite 104ca Bragg peak is also random, but fluctuates strongly as function of . The data thus suggests that there are many quartz crystallites, but only a few calcite crystallites. For both c) and d), the scattering at  ~ 45º is missing due to the beam stop support.

Meridional and equatorial descriptors
The idea coined "meridional descriptor" is essentially to integrate over a certain region of Qspace in such a way that the scattered intensity is minimally affected under sample rotation. Fig. S2a illustrates the Ewald sphere when considering diffraction from a single set of planes, with Ghkl and -Ghkl denoting a generic pair of reciprocal lattice vectors. In the general case when the sample is rotated, Ghkl will intersect the Ewald sphere for some sample rotation angles ω, causing diffraction peaks that can be observed on the detector. It is clear that Bragg peaks with preferred orientation in the x-z plane will change rapidly in intensity as the sample rotates. However, if Ghkl and -Ghkl point essentially in the meridional ("vertical"; y) direction, parallel to the tomography axis, the diffracted intensity will to a first approximation be independent of ω. Because the photon energy used is high, the Ewald sphere is almost flat, and the Bragg angles θ small. If θ is comparable or smaller than the width of the orientation distribution Δ, the vertical diffraction would be expected to be essentially constant when the sample is rotated, and thus we can define the meridional descriptor I|| as mer hkl || hkl measured , , , Formally, I|| differs from Iisotropic by integrating  over only the meridional sectors of the detector. For example, for the 002HA peak at θ = 1.2º for the energy used in this experiment, I() has a width Δ  50º, which is much larger than the corresponding Bragg angle, implying that if G002 is vertical, it will give essentially the same scattered intensity on the detector regardless of the sample orientation . The meridional descriptor proved essential for retrieving the orientation-dependent tomograms. It should be noted that the method described in this article, using the meridional descriptor, only provides estimates of the vectorial tomograms. In order to obtain a vectorial tomogram in a mathematically correct way, it is necessary to measure all three components of a lattice vector Ghkl for each voxel. The meridional descriptor does not measure these three components, but returns a measure of the statistical distribution of crystallites orientation within each voxel. This further implies that measuring samples with three orthogonal orientations is not sufficient for true vectorial tomography, but can provide useful qualitative results.
Similarly, the equatorial descriptor I⊥ is defined for scattering to the equatorial detector sectors as eq hkl hkl measured ,, For crystallites with a preferred orientation, I⊥ from a sample volume element will generally change during sample rotation, and is therefore in general not suited for tomographic reconstruction using the FBP algorithm. However, when assumed that the crystallite orientation distributions are wide, and that the chosen Bragg peak for reconstruction is aligned along the tomography axis, one can obtain information about crystallite orientation, as shown for the samples from Discosauriscus austriacus and Eusthenopteron foordi.
As shown in Fig. S2b-d, the different minerals constituting the sample were found to exhibit qualitatively different orientation distributions. Nevertheless, 3D compositional tomograms for all minerals could be generated with few artefacts by using the isotropic "powder" intensity descriptor (Equation S1 ). The 3D fossil bone composition is shown with 3D XRD-CT tomograms in Fig. 4. It is particularly satisfying to note how the different minerals constituting the sample fill complementary spatial regions. Unlike the other compounds, which form rather connected regions extending to the sample surface, the barite crystallites were found to be distributed with isolated regions solely in the sample interior.

Minerals present
Diffraction patterns representative of large regions of the sample were obtained by summing the diffraction patterns for many different sample positions and orientations, cf. Fig. S3 for the azimuthally averaged I(Q). Figure S3 A 1D intensity global pattern for the humerus of Eusthenopteron foordi obtained by summing many individual azimuthally averaged diffraction patterns, used for identifying Bragg peaks suitable as signal for CT reconstruction. The most intense Bragg peaks belonging to the different minerals have been indicated with coloured lines. The background is seen to be essentially constant.
Bragg peak positions expected for the different minerals contained in the sample have been indicated with coloured lines in Fig. S3. The Bragg peak positions observed from the experimental data were compared with tabulated data [4][5][6][7][8] , and all major Bragg peaks could be identified as originating from the listed five minerals. The relatively large pixel size of the detector combined with the small scattering angles caused a significant number of overlapping peaks, and only a few of the observed Bragg peaks could be extracted reliably and isolated for tomographic reconstruction, as listed in Table S1. The minerals in the sample, quartz, calcite barite and pyrite, were seen to exhibit qualitative differences in the scattering patterns in terms of preferred orientation, size and quantity. Large crystallites necessarily cause streak artefacts in the XRD-CT tomograms, but it was possible to reduce these reconstruction artefacts by combining diffraction intensities from several Bragg peaks, as well as by combining the results from the three different datasets measured. Note that for pyrite, all except one observed Bragg peak were overlapping with peaks from HA and barite, which made it difficult to provide reliable pyrite tomograms. For the HA crystallites, the reconstructed tomograms were quantitatively similar regardless of which intensity descriptor that was used as input to the FBP. As shown with an example in Fig. S2b, the scattering signal from HA was smoothly varying with sample position and orientation, consistent with the presence of numerous crystallites, presumably of small size compared to the sample volume. A significant fraction of the crystallites was isotropically oriented, and using the isotropic descriptor resulted in tomograms with few streak artefacts.
For quartz, Fig. S2c shows that the crystallites diffract approximately equally to all  for the 100qu Bragg peak, suggesting that the quartz crystallites are isotropically oriented and that diffraction from a large number of crystallites is observed. The observed intensity from calcite fluctuates strongly with , as shown in Fig. S2d for the 202ca Bragg peak. As for quartz, the calcite crystallites appear randomly oriented, but the Bragg intensities fluctuate strongly with , which indicates the presence of relatively few crystallites. The distribution of observed intensity from the barite 311ba Bragg peak (not shown) exhibits sharp and strong peaks scattering randomly distributed in , also suggesting the presence of rather few crystallites.