# Automation of the Peak Fitting Method in Bone FTIR Microspectroscopy Spectrum Analysis: Human and Mice Bone Study

Authors: Marc Gardegaront; Delphine Farlay; Olivier Peyruchaud; Hélène Follet
Journal: Journal of Spectroscopy
 (2018)

---

## Abstract

FTIR microspectroscopy (FTIRM) is a commonly used nondestructive method to characterise thin bone sections. However, spectrum analysis methods are often highly sensitive to small variations (e.g., boundary limits), thus implying a time-consuming and redundant analysis process. To solve this issue, software has been developed based on several algorithms to automate the analysis. Furthermore, a rigorous framework has been established concerning the peak fitting method to obtain the systematic best potential solution. Validation of the automatic method has been performed by comparison with the manual method. Results and validation proved the reliability of the automatic process. The developed algorithms provide the means necessary to fully compare the results between bone FTIRM studies and between different laboratories.


---

## Body

## 1. Introduction

Bone analysis methods are quite varied [1–3], ranging from the noninvasive, e.g., X-ray scanning, to the destructive, e.g., mechanical rupture tests. FTIR microspectroscopy (FTIRM) is a promising field which characterises bone samples, as it is able to provide information on the chemical composition and structural organisation of the analysed matter without damaging it on a local level of a thin section [4–9]. This method combines optic, quantum mechanics, chemistry, and biology knowledge, which leads to a complex understanding that might have an impact on the conclusion of an FTIR study [10]. In addition, the spectrum analysis is laborious if done “manually.”Several different methods are currently used to analyse FTIRM spectra (peak fitting method [11–13], second derivative method [13, 14], etc.), and each has its advantages and limitations; moreover, none are considered to be a gold standard. We chose to stick with the peak fitting method used in [11, 12]. The second derivative method, although quite promising due to its objective and purely mathematical analysis, is too sensitive to noise (producing peaks that are not there initially) and quite inadequate in the case of peaks that are too close and with different full width at half maximum (FWHM) (the wider peak is then masked by the thinner one). This is why the second derivative method is not chosen (for more explanation, cf. Appendix C).Traditionally, the peak fitting of a bone spectrum is done by selecting a region of interest (ROI) in the bone spectrum, followed by fitting the curve with a sum of Gaussians. Previously, this fitting was done by assuming both position and characteristics of the peaks (height and FWHM), giving a starting state for the least-squares fitting method. The result of the fitting process was then left to the operator for visual verification. As one can easily imagine, this technique proved to be fastidious, as correct results were typically obtained through several tries. Indeed, the regression may sometimes be highly sensitive and be executed multiple times for one ROI in order to find the correct fitting result.In a bone spectrum (Figure1), six distinct regions are identified as the vibrational response of different constituents:(i)
The amide (I, II, and III) areas (1300, 1700) cm−1 corresponding to the response of the organic matrix in the bone tissue (with also a ν3CO3 contribution at around 1400 cm−1)(ii)
Theν1ν3PO4 area (900, 1200) cm−1 corresponding to the symmetric and antisymmetric stretching of the phosphates(iii)
Theν2CO3 area (800, 900) cm−1 corresponding to the symmetrical stretching of the carbonates (CO3), which also contains a large HPO4 contribution, which is removed before analysis by using a wisely placed baseline correction(iv)
Theν4PO4 area (500, 650) cm−1 corresponding to the antisymmetric deformation of the phosphates (PO4).Figure 1
Typical FTIRM spectrum of human bone between the 450 cm−1 and 4000 cm−1 boundary. Amides I, II, and III represent the organic matrix of the bone, while the ν1ν3PO4, ν4PO4, and ν2CO3 are part of the mineral matrix.We automated the methods used in [9, 11, 14, 15], which are biologically relevant methods (not meant to be spectroscopically accurate), to obtain the five parameters within these regions to characterise the sample:Directly measured on the spectra are as follows:(i)
Mineral-to-organic ratio (mineral/matrix): ratio of theν1ν3PO4 area (910, 1184) cm−1 over the amide I area (1592, 1730) cm−1 [14](ii)
Carbonation: ratio of theν2CO3 area (862, 894) cm−1 over the ν1ν3PO4 area (910, 1184) cm−1 [14].Postpeak fitting is given as follows(i)
Crystallinity Index: inverse of the FWHM of the ∼604 cm−1 peak [15](ii)
Mineral Maturity (MM): ratio of the apatitic (∼1030 cm−1 peak) over nonapatitic (∼1110 cm−1 peak) part of the ν1ν3PO4 area [15](iii)
Collagen Maturity (Col Mat): ∼1660 cm−1 peak over ∼1690 cm−1 peak (part of the amide I area) [5, 11].The manual peak fitting method might have induced subjectivity and imprecision, as there were no strict means to analyse bone spectra. Also, these methods are not always spectroscopically accurate (i.e., there are not always the right number of peaks corresponding to the exact composition of the studied material); however, these methods have been validated and are sufficient to conclude the proportion of compounds in the analysed bone.These methods might be improved by the careful study and translation of the analysis part into an algorithm implemented into custom software.The aim of this technical note was to offer a solution to the automation process of the peak fitting method used in FTIR spectrum analysis of the bone. Software, i.e., LYOS spectrum analysis, written in Python, has been developed to fully automate the analysis process of healthy human and thin mouse bone section’s infrared spectra, even though the software can be used to analyse other types of spectrum.


## 2. Materials and Methods

The validation process has been done using the following:(1)
Fifty-one trabecular thin sections of human L2 vertebra from 51 donors (French body donation to science program, declaration number: DC-2015-2357; Laboratory of Anatomy, Faculty of Medicine Lyon Est, University of Lyon, France) (20 spectra per sample at 30 × 100µm2 spatial resolution)(2)
Four trabecular and cortical thin sections of human iliac crest from four donors (Laboratoire d’anatomie, Université Claude Bernard, Lyon, France) (78 spectra at 50 × 50µm2 spatial resolution)(3)
Twelve cortical thin sections of human femur from 12 subjects (MULTIPS: ANR-13-BS09-0006-01) (12 spectra at 150 × 150µm2 spatial resolution)(4)
One trabecular and cortical thin section of mouse tibia from one mouse (Balb/C E2A-55 #DR2015-58 MENESR: #APAFIS# 39132015121 0903612v2, Lysbone: ANR-15- CE14-0010) (10 spectra at 30 × 30µm2 spatial resolution).After dissection, samples were fixed in 70% ethanol for two weeks, dehydrated for 48 h in 100% ethanol, substituted for 48 h in methylcyclohexane (solvent that removes lipids in the samples), and then embedded in methyl methacrylate (MMA). The sections of 2µm thickness were obtained with a microtome Polycut E (Reichert-Jung, Leica, Germany) and stored between two slides at room temperature.The raw spectra were collected at 1 cm−1 resolution and averaged by 40 scans per scanned region in transmission mode with a Perkin–Elmer GXII AutoIMAGE microscope (Norwalk, CT, USA) equipped with a wide-band detector (mercury-cadmium-telluride) (7800–400 cm−1). The instrument used a Cassegrain objective of numerical aperture 0.6. The system has a spatial resolution of 10 µm at typical midinfrared wavelengths. Contribution from air was subtracted from the original spectrum.The commercial software Spectrum (v6.2.0.0055) provided by Perkin–Elmer and Grams/AI 8.0 Spectroscopy Software provided by the Thermo Electron Corporation offer various preprocessing options and customisation for the automatic peak fitting process. The mandatory methods used in these software have been implemented and improved for automation in the custom Python software and are described below, such as a transformation from transmittance to absorbance values, a spectrum component removal (to extract the components of interest from a raw spectrum), and a baseline correction method (to remove the high-scale effects of a bone FTIRM measure).


### 2.1. Transmittance to Absorbance

The FTIRM spectrum of a thin section of the bone is obtained in transmittance. However, the spectrum analysis is conventionally done on an absorbance spectrum. The transmittance spectrum (T, in percent) is then transformed into its absorbance (A) form by the following function:(1)A=−log10T100.


### 2.2. Spectrum Component Removal

Thin bone sections could not be realised without an embedding medium (methyl methacrylate, MMA). This medium creates a potential infrared-reactive component in the spectrum that must be removed in order to obtain the bone-only spectrum. This component presents a high absorption peak around 1730 cm−1 (Figure 2). The first step is to scan the MMA to have its pure spectrum in the same range as the bone spectrum (450, 4000) cm−1. Then, a custom algorithm of optimisation subtracts the MMA spectrum from the bone spectrum (cf. Appendix A). This process is done by selecting the boundaries of the highest MMA peak (1826, 1706) cm−1 between which is computed with the second derivative of the bone spectrum (corresponding to the curvature of the spectrum). An iterative loop consecutively guesses a coefficient by which the whole MMA spectrum is multiplied, and the subtraction result is evaluated on the convexity of the shape of the spectrum. The convexity of a curve is mathematically defined by the positivity of the second derivative. The best result is the one that has the minimum negative part. Each loop will check the negative part of the second derivative of the bone-embedded spectrum and try to minimise it by subtracting the weighted MMA-only spectrum (Figure 3). The algorithm can be summed up as follows:Findyi which ∫fyi is minimumwithyi: the spectrum subtracted by the ith weighted MMA-only spectrumandfyi=0,ifyi>0,yi,ifyi≤0.Figure 2
Typical sample of bone embedded in MMA (left) with its respective spectrum (top right) and the pure MMA spectrum (bottom right). A: cortical bone (human femur) embedded in MMA. B: methyl methacrylate (MMA). The MMA spectrum is taken outside the bone region of the sample.Figure 3
Visual representation of the algorithm. Evolution of the negative part of the second derivative of the bone spectrum between the MMA peak of interest boundaries (1706 cm−1 and 1826 cm−1) through the iterations. Each iteration corresponds to the subtraction of the MMA-embedded bone spectrum minus the weighted MMA spectrum (the weight changes with each iteration). Note. An iteration is not always better than the previous one (e.g., iterations 1 and 2), but it tends nonetheless to converge to the expected result (e.g., iteration 22), which is the convexity of the spectrum between the 1706 cm−1 and 1826 cm−1 boundaries. The final result of the MMA removal can be observed in Figure 1 (the MMA contribution, including the 1730 cm−1 peak, in the bone spectrum is removed).


### 2.3. Baseline Correction

The baseline of the bone in the (450, 4000) cm−1 is convex. The purpose is to realign this line horizontally. A multilinear baseline is chosen as an approximation of the effects of the interferences (beam spreading, sample heating, etc). Based on the convexity hypothesis, the algorithm used to find the baseline will first detect the local minima using a Savitzky–Golay smoothing [16, 17] (keeping intensity and position of the local extrema of the raw spectrum) and draw the segments between two minima containing the most minima above it while forbidding the presence of minima beneath it (Figure 4).Figure 4
Graphical representation of the baseline creation. Red dotes are local minima, and the black line is the created baseline.


### 2.4. Peak Fitting

The peak fitting method is one of the methods used in the conventional FTIRM spectrum analysis [18] (least-squares method with trust region reflective algorithm to find the best suited parameters). The curve fitting function used is curve_fit from the scipy.optimize module version 1.0.0. Usually, only a first guess state for each Gaussian peak is given for the manual method, allowing the peaks to place everywhere with any characteristics (position, FWHM, and height). This first guess-only method often induces inaccurate results and forces the operator to execute the method multiple times in order to find a valid fit. The automation process is implemented by adding boundaries to find each peak’s parameter. These boundaries are carefully chosen for the solution to be unique (restricting the plurality of solutions induced by the manual method with first state guessing only). Ideally, boundaries should be chosen considering the characteristics of each sample (temperature, internal constraints, crystal properties, etc.). However, which characteristic and how each influences a peak’s parameter cannot be answered precisely with our current knowledge. Thus, a more empirical approach has been chosen to determine the boundaries of each peak. Based on a set of more than a thousand spectra considered to be well peak fitted, some schemes have been observed. Eight peaks were found in the amide area, five peaks in the ν1ν3PO4 area, and four peaks in the ν4PO4 area. Boundaries for each peak’s parameters have been established and are registered in Table 1. For each area, the first guess of each peak located at a specific wave number is set at 90% of the absorbance in height and at 5 cm−1 in FWHM. The ranges for the peak fitting of each area are chosen to give mathematically correct results. They are related to the biologically relevant range given in the introduction part but generally need to be wider, as mathematical information is present beyond the biologically relevant range (for example, the amide area range is (1300, 1700) cm−1 and the mathematical range used is (1300, 1787) cm−1, the same for the ν1ν3PO4 and ν4PO4). The algorithm can be summed up as follows:Find thejth Gaussians parameters Gj,x0j,ωjfor whichfx minimizes ∑iyi−fxi2withfx=∑jGj∗e−x−x0j/ωj2andGj: amplitude of the jth Gaussian,x0j: position of the jth Gaussian,ωj: FWHM/22 of the jth GaussianTable 1
Boundaries of each peak used for the peak fitting method. Theν4PO4 area contains 4 peaks, the ν1ν3PO4 area contains 5 peaks, and the amide area contains 8 peaks. Each peak’s parameter (height, FWHM, and position) is searched between the defined minimum and maximum boundaries, expressed in the table, by the fitting algorithm.
ν
4PO4 area (531, 659) cm−1
ν
1
ν
3PO4 area (910, 1184) cm−1
Amide area (1300, 1787) cm−1
First guess
Position
540
565
575
604
962
1030
1060
1070
1110
1340
1420
1450
1550
1600
1633
1660
1690
Height
Minimum
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Maximum
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
FWHM (cm−1)
Minimum
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Maximum
60
50
40
∞
70
70
40
40
70
∞
∞
∞
∞
30
∞
∞
∞
Position (cm−1)
Minimum
530
555
565
−∞
932
1000
1045
1055
1080
1310
1390
1420
1520
1585
1633
1631
1690
Maximum
550
575
585
+∞
992
1060
1075
1085
1140
1370
1450
1480
1580
1615
1633
1689
1690
∞: infinity.with respect to the constraints to each Gaussian’s parameters explicited in Table1.


#### 2.4.1. The Amide Area (1300, 1787) cm−1

The area of interest inside the amide part of the spectrum lies around the (1600, 1730) cm−1 boundary (Figure 5). Even if (1300, 1600) cm−1 is not used in the measured parameters, it was kept to accurately fit the (1300, 1730) cm−1 area. The main peaks are the dashed-blue peak (∼1690 cm−1) and the dashed-purple peak (∼1660 cm−1). The 1660 cm−1 peak is highly influenced by the 1633 cm−1 (dashed orange), which is as well influenced by the 1600 cm−1 peak (full green peak). In bone, the 1633 cm−1 peak must be higher than the local minima at 1600 cm−1, and the height of the 1600 cm−1 peak greatly influences the height of the 1633 cm−1. All the parameters are constrained with static values (unchanged regarding the spectrum) (cf. Table 1). It may happen that the MMA peak is not correctly removed due to small differences of the shape of the MMA peak between the pure MMA spectrum and the MMA component in the bone + MMA spectrum. In order to minimise the implications of the residual MMA peak, an additional peak is included for the peak fitting method, which is located at 1745 cm−1. This is an automatic process realised when the algorithm detects a significant peak between 1730 cm−1 and 1760 cm−1 (i.e., when the positive area of the derivate between those values is superior to the half of the negative area of the derivate).Figure 5
Original spectrum in the amide area (dark blue line) recovered by the Gaussian sum fitting (dotted-black line) and the 8 individual peaks beneath (dashed and full lines) + 1 peak for the MMA (1730 cm−1). The two peaks directly used to compute the Collagen Maturity parameter are the 1660 cm−1 peak (dashed-purple line) and the 1690 cm−1 peak (dashed-blue one).


#### 2.4.2. Theν1ν3PO4 Area (910, 1184) cm−1

The two peaks of interest here are the 1030 and 1110 cm−1 peak (full purple and green) (Figure 6). These two peaks are mainly influenced by the position of the peaks between 1060 and 1070 cm−1. Every parameter for this area determination is a static value (cf. Table 1).Figure 6
Original spectrum in theν1ν3PO4 area (dark blue line) recovered by the Gaussian sum fitting (dotted-black line) and the 5 individual peaks beneath (full lines). The two peaks directly used to compute the Mineral Maturity parameter are the 1030 cm−1 peak (purple line) and the 1110 cm−1 peak (dashed-green one).


#### 2.4.3. Theν4PO4 Area (531, 659) cm−1

Only the 604 cm−1 peak (full yellow) is used to characterise the bone in the ν4PO4 area (Figure 7). This peak is mainly influenced by the 575 cm−1 peak. In bone, the 575 cm−1 peak’s height must be inferior to the 565 cm−1 peak to avoid the 575 cm−1 peak to spread too much on the 604 cm−1 peak area. Every parameter for this area is a static value (cf. Table 1).Figure 7
Original spectrum in theν4PO4 area (dark blue line) recovered by the Gaussian sum fitting (dotted-black line) and the 4 individual peaks beneath (full lines). The peak used to compute the Crystallinity parameter is 604 cm−1 (yellow line).


#### 2.4.4. Statistical Analysis

The validation of the spectrum analysis methods has been done using the root mean square of the coefficient of variation to assess the variation in manual methods, with the Bland–Altman method to compare manual and automatic results.


## 2.1. Transmittance to Absorbance

The FTIRM spectrum of a thin section of the bone is obtained in transmittance. However, the spectrum analysis is conventionally done on an absorbance spectrum. The transmittance spectrum (T, in percent) is then transformed into its absorbance (A) form by the following function:(1)A=−log10T100.


## 2.2. Spectrum Component Removal

Thin bone sections could not be realised without an embedding medium (methyl methacrylate, MMA). This medium creates a potential infrared-reactive component in the spectrum that must be removed in order to obtain the bone-only spectrum. This component presents a high absorption peak around 1730 cm−1 (Figure 2). The first step is to scan the MMA to have its pure spectrum in the same range as the bone spectrum (450, 4000) cm−1. Then, a custom algorithm of optimisation subtracts the MMA spectrum from the bone spectrum (cf. Appendix A). This process is done by selecting the boundaries of the highest MMA peak (1826, 1706) cm−1 between which is computed with the second derivative of the bone spectrum (corresponding to the curvature of the spectrum). An iterative loop consecutively guesses a coefficient by which the whole MMA spectrum is multiplied, and the subtraction result is evaluated on the convexity of the shape of the spectrum. The convexity of a curve is mathematically defined by the positivity of the second derivative. The best result is the one that has the minimum negative part. Each loop will check the negative part of the second derivative of the bone-embedded spectrum and try to minimise it by subtracting the weighted MMA-only spectrum (Figure 3). The algorithm can be summed up as follows:Findyi which ∫fyi is minimumwithyi: the spectrum subtracted by the ith weighted MMA-only spectrumandfyi=0,ifyi>0,yi,ifyi≤0.Figure 2
Typical sample of bone embedded in MMA (left) with its respective spectrum (top right) and the pure MMA spectrum (bottom right). A: cortical bone (human femur) embedded in MMA. B: methyl methacrylate (MMA). The MMA spectrum is taken outside the bone region of the sample.Figure 3
Visual representation of the algorithm. Evolution of the negative part of the second derivative of the bone spectrum between the MMA peak of interest boundaries (1706 cm−1 and 1826 cm−1) through the iterations. Each iteration corresponds to the subtraction of the MMA-embedded bone spectrum minus the weighted MMA spectrum (the weight changes with each iteration). Note. An iteration is not always better than the previous one (e.g., iterations 1 and 2), but it tends nonetheless to converge to the expected result (e.g., iteration 22), which is the convexity of the spectrum between the 1706 cm−1 and 1826 cm−1 boundaries. The final result of the MMA removal can be observed in Figure 1 (the MMA contribution, including the 1730 cm−1 peak, in the bone spectrum is removed).


## 2.3. Baseline Correction

The baseline of the bone in the (450, 4000) cm−1 is convex. The purpose is to realign this line horizontally. A multilinear baseline is chosen as an approximation of the effects of the interferences (beam spreading, sample heating, etc). Based on the convexity hypothesis, the algorithm used to find the baseline will first detect the local minima using a Savitzky–Golay smoothing [16, 17] (keeping intensity and position of the local extrema of the raw spectrum) and draw the segments between two minima containing the most minima above it while forbidding the presence of minima beneath it (Figure 4).Figure 4
Graphical representation of the baseline creation. Red dotes are local minima, and the black line is the created baseline.


## 2.4. Peak Fitting

The peak fitting method is one of the methods used in the conventional FTIRM spectrum analysis [18] (least-squares method with trust region reflective algorithm to find the best suited parameters). The curve fitting function used is curve_fit from the scipy.optimize module version 1.0.0. Usually, only a first guess state for each Gaussian peak is given for the manual method, allowing the peaks to place everywhere with any characteristics (position, FWHM, and height). This first guess-only method often induces inaccurate results and forces the operator to execute the method multiple times in order to find a valid fit. The automation process is implemented by adding boundaries to find each peak’s parameter. These boundaries are carefully chosen for the solution to be unique (restricting the plurality of solutions induced by the manual method with first state guessing only). Ideally, boundaries should be chosen considering the characteristics of each sample (temperature, internal constraints, crystal properties, etc.). However, which characteristic and how each influences a peak’s parameter cannot be answered precisely with our current knowledge. Thus, a more empirical approach has been chosen to determine the boundaries of each peak. Based on a set of more than a thousand spectra considered to be well peak fitted, some schemes have been observed. Eight peaks were found in the amide area, five peaks in the ν1ν3PO4 area, and four peaks in the ν4PO4 area. Boundaries for each peak’s parameters have been established and are registered in Table 1. For each area, the first guess of each peak located at a specific wave number is set at 90% of the absorbance in height and at 5 cm−1 in FWHM. The ranges for the peak fitting of each area are chosen to give mathematically correct results. They are related to the biologically relevant range given in the introduction part but generally need to be wider, as mathematical information is present beyond the biologically relevant range (for example, the amide area range is (1300, 1700) cm−1 and the mathematical range used is (1300, 1787) cm−1, the same for the ν1ν3PO4 and ν4PO4). The algorithm can be summed up as follows:Find thejth Gaussians parameters Gj,x0j,ωjfor whichfx minimizes ∑iyi−fxi2withfx=∑jGj∗e−x−x0j/ωj2andGj: amplitude of the jth Gaussian,x0j: position of the jth Gaussian,ωj: FWHM/22 of the jth GaussianTable 1
Boundaries of each peak used for the peak fitting method. Theν4PO4 area contains 4 peaks, the ν1ν3PO4 area contains 5 peaks, and the amide area contains 8 peaks. Each peak’s parameter (height, FWHM, and position) is searched between the defined minimum and maximum boundaries, expressed in the table, by the fitting algorithm.
ν
4PO4 area (531, 659) cm−1
ν
1
ν
3PO4 area (910, 1184) cm−1
Amide area (1300, 1787) cm−1
First guess
Position
540
565
575
604
962
1030
1060
1070
1110
1340
1420
1450
1550
1600
1633
1660
1690
Height
Minimum
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Maximum
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
∞
FWHM (cm−1)
Minimum
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Maximum
60
50
40
∞
70
70
40
40
70
∞
∞
∞
∞
30
∞
∞
∞
Position (cm−1)
Minimum
530
555
565
−∞
932
1000
1045
1055
1080
1310
1390
1420
1520
1585
1633
1631
1690
Maximum
550
575
585
+∞
992
1060
1075
1085
1140
1370
1450
1480
1580
1615
1633
1689
1690
∞: infinity.with respect to the constraints to each Gaussian’s parameters explicited in Table1.


### 2.4.1. The Amide Area (1300, 1787) cm−1

The area of interest inside the amide part of the spectrum lies around the (1600, 1730) cm−1 boundary (Figure 5). Even if (1300, 1600) cm−1 is not used in the measured parameters, it was kept to accurately fit the (1300, 1730) cm−1 area. The main peaks are the dashed-blue peak (∼1690 cm−1) and the dashed-purple peak (∼1660 cm−1). The 1660 cm−1 peak is highly influenced by the 1633 cm−1 (dashed orange), which is as well influenced by the 1600 cm−1 peak (full green peak). In bone, the 1633 cm−1 peak must be higher than the local minima at 1600 cm−1, and the height of the 1600 cm−1 peak greatly influences the height of the 1633 cm−1. All the parameters are constrained with static values (unchanged regarding the spectrum) (cf. Table 1). It may happen that the MMA peak is not correctly removed due to small differences of the shape of the MMA peak between the pure MMA spectrum and the MMA component in the bone + MMA spectrum. In order to minimise the implications of the residual MMA peak, an additional peak is included for the peak fitting method, which is located at 1745 cm−1. This is an automatic process realised when the algorithm detects a significant peak between 1730 cm−1 and 1760 cm−1 (i.e., when the positive area of the derivate between those values is superior to the half of the negative area of the derivate).Figure 5
Original spectrum in the amide area (dark blue line) recovered by the Gaussian sum fitting (dotted-black line) and the 8 individual peaks beneath (dashed and full lines) + 1 peak for the MMA (1730 cm−1). The two peaks directly used to compute the Collagen Maturity parameter are the 1660 cm−1 peak (dashed-purple line) and the 1690 cm−1 peak (dashed-blue one).


### 2.4.2. Theν1ν3PO4 Area (910, 1184) cm−1

The two peaks of interest here are the 1030 and 1110 cm−1 peak (full purple and green) (Figure 6). These two peaks are mainly influenced by the position of the peaks between 1060 and 1070 cm−1. Every parameter for this area determination is a static value (cf. Table 1).Figure 6
Original spectrum in theν1ν3PO4 area (dark blue line) recovered by the Gaussian sum fitting (dotted-black line) and the 5 individual peaks beneath (full lines). The two peaks directly used to compute the Mineral Maturity parameter are the 1030 cm−1 peak (purple line) and the 1110 cm−1 peak (dashed-green one).


### 2.4.3. Theν4PO4 Area (531, 659) cm−1

Only the 604 cm−1 peak (full yellow) is used to characterise the bone in the ν4PO4 area (Figure 7). This peak is mainly influenced by the 575 cm−1 peak. In bone, the 575 cm−1 peak’s height must be inferior to the 565 cm−1 peak to avoid the 575 cm−1 peak to spread too much on the 604 cm−1 peak area. Every parameter for this area is a static value (cf. Table 1).Figure 7
Original spectrum in theν4PO4 area (dark blue line) recovered by the Gaussian sum fitting (dotted-black line) and the 4 individual peaks beneath (full lines). The peak used to compute the Crystallinity parameter is 604 cm−1 (yellow line).


### 2.4.4. Statistical Analysis

The validation of the spectrum analysis methods has been done using the root mean square of the coefficient of variation to assess the variation in manual methods, with the Bland–Altman method to compare manual and automatic results.


## 2.4.1. The Amide Area (1300, 1787) cm−1

The area of interest inside the amide part of the spectrum lies around the (1600, 1730) cm−1 boundary (Figure 5). Even if (1300, 1600) cm−1 is not used in the measured parameters, it was kept to accurately fit the (1300, 1730) cm−1 area. The main peaks are the dashed-blue peak (∼1690 cm−1) and the dashed-purple peak (∼1660 cm−1). The 1660 cm−1 peak is highly influenced by the 1633 cm−1 (dashed orange), which is as well influenced by the 1600 cm−1 peak (full green peak). In bone, the 1633 cm−1 peak must be higher than the local minima at 1600 cm−1, and the height of the 1600 cm−1 peak greatly influences the height of the 1633 cm−1. All the parameters are constrained with static values (unchanged regarding the spectrum) (cf. Table 1). It may happen that the MMA peak is not correctly removed due to small differences of the shape of the MMA peak between the pure MMA spectrum and the MMA component in the bone + MMA spectrum. In order to minimise the implications of the residual MMA peak, an additional peak is included for the peak fitting method, which is located at 1745 cm−1. This is an automatic process realised when the algorithm detects a significant peak between 1730 cm−1 and 1760 cm−1 (i.e., when the positive area of the derivate between those values is superior to the half of the negative area of the derivate).Figure 5
Original spectrum in the amide area (dark blue line) recovered by the Gaussian sum fitting (dotted-black line) and the 8 individual peaks beneath (dashed and full lines) + 1 peak for the MMA (1730 cm−1). The two peaks directly used to compute the Collagen Maturity parameter are the 1660 cm−1 peak (dashed-purple line) and the 1690 cm−1 peak (dashed-blue one).


## 2.4.2. Theν1ν3PO4 Area (910, 1184) cm−1

The two peaks of interest here are the 1030 and 1110 cm−1 peak (full purple and green) (Figure 6). These two peaks are mainly influenced by the position of the peaks between 1060 and 1070 cm−1. Every parameter for this area determination is a static value (cf. Table 1).Figure 6
Original spectrum in theν1ν3PO4 area (dark blue line) recovered by the Gaussian sum fitting (dotted-black line) and the 5 individual peaks beneath (full lines). The two peaks directly used to compute the Mineral Maturity parameter are the 1030 cm−1 peak (purple line) and the 1110 cm−1 peak (dashed-green one).


## 2.4.3. Theν4PO4 Area (531, 659) cm−1

Only the 604 cm−1 peak (full yellow) is used to characterise the bone in the ν4PO4 area (Figure 7). This peak is mainly influenced by the 575 cm−1 peak. In bone, the 575 cm−1 peak’s height must be inferior to the 565 cm−1 peak to avoid the 575 cm−1 peak to spread too much on the 604 cm−1 peak area. Every parameter for this area is a static value (cf. Table 1).Figure 7
Original spectrum in theν4PO4 area (dark blue line) recovered by the Gaussian sum fitting (dotted-black line) and the 4 individual peaks beneath (full lines). The peak used to compute the Crystallinity parameter is 604 cm−1 (yellow line).


## 2.4.4. Statistical Analysis

The validation of the spectrum analysis methods has been done using the root mean square of the coefficient of variation to assess the variation in manual methods, with the Bland–Altman method to compare manual and automatic results.


## 3. Results and Discussion

### 3.1. Validation

The first step of the validation process is done by analysing the precision of the manual peak fitting method. Two types of precision are measured: (1) the intraoperator precision that measures the ability of one operator to obtain the same result multiple times on the same spectrum and (2) the interoperator precision that measures the ability to get the same result from two different operators on the same spectrum.The first operator (operator 1) is experienced with the bone FTIRM peak fitting, and is considered the reference. The second (operator 2) is an inexperienced operator who has been given the instructions to execute the peak fitting analysis.Five different bone spectra have been peak fitted with the manual method two times by each operator (cf. Table2). For the manual peak fitting, the root mean square of the coefficient of variation shows that results are operator-dependant.Table 2
Differences of the manual peak fitting analysis intraoperator, interoperator, maximum absolute error, and variation of results using the automatic method for each parameter, based on the analysis of 5 bone spectra. Operator 1 is an experienced operator, and operator 2 is an inexperienced operator to whom has been given the instructions to analyse the spectra.
Crystallinity Index
Mineral Maturity
Collagen Maturity
RMS of coefficient of variation (operator 1) (%)
0.7
1.4
0.7
RMS of coefficient of variation (operator 2) (%)
1.7
0.4
1.0
RMS of coefficient of variation (interoperator) (%)
2.0
1.2
0.7
Maximum absolute error
0.0013 cm
0.07
0.11
Variation of results using automatic method
No variation
No variation
No variation
RMS = root mean square.Then, we compared both automatic and manual methods based on the spectra extracted from the human samples and mouse sample. The Bland–Altman graph plots the means and the differences between automatic and manual methods measured by a human operator. For the human bones, Figure8 shows, at first sight, an almost null shift for every parameter. The ratio of valid results (i.e., results between the mean ± 1.96 ∗ SD) is 92% for the Crystallinity Index, 93% for the Mineral Maturity, and 92% for the Collagen Maturity.Figure 8
The Bland–Altman plot of the human bone spectra (manual method-automatic method values). The straight line is the mean of the differences. The dashed line is the mean of the differences plus or minus 1.96 ∗ SD, representing the limit of agreement. 92% of Crystallinity Index values, 93% of Mineral Maturity values, and 92% of Collagen Maturity values are within the dashed lines. Shifts of 0.0007 cm for the Crystallinity Index, 0.11 for the Mineral Maturity, and 0.17 for the Collagen Maturity are observed.
(a)
(b)
(c)For the mice tibia (Figure9), the shift is also almost null for every parameter. The ratio of valid results under the aforementioned hypothesis is 90% of valid results for the Crystallinity Index, 100% for the Mineral Maturity, and 100% for the Collagen Maturity.Figure 9
The Bland–Altman plot of the mice bone spectra (manual method-automatic method values). The straight line is the mean of the differences. The dashed line is the mean of the differences plus or minus 1.96 ∗ SD, representing the limit of agreement. 90% of Crystallinity Index values, 100% of Mineral Maturity values, and 100% of Collagen Maturity values are within the dashed lines. Shifts of −0.0009 cm for the Crystallinity Index, −0.01 for the Mineral Maturity, and 0.21 for the Collagen Maturity are observed.
(a)
(b)
(c)


### 3.2. Discussion

Gaussians have been long used to establish methods and peaks of interests in human bone, meaning changing the fitting function would end in a change of the peaks of interest. Thus, to maintain consistency with previous studies, Gaussians are kept for the fitting method. However, it should be noted that Lorentzians or Voigt functions would be more accurate regarding the quantum reality of the spectra [19].The development of the automatic method not only allowed a significant gain of time but also suppressed the interoperator and intraoperator variability.Currently, there is no gold standard to check the exactitude of the results obtained. After careful examination of the automatic peak fitting of these spectra, it has been found that the results obtained with the automatic method are valid compared with the analysis done by the expert using the manual method.The observed shifts are mainly due to the choice of boundaries within which the peak fitting method is applied. Indeed, there are no strict boundaries required for the manual peak fitting of theν4PO4 area. This induces the variation of the shift, as shown in Appendix B. This problem is solved by the automatic method that fixes those boundaries, chosen after careful analysis of several sets of boundaries.We could argue that it is not the most exact method in terms of localisation of the peaks, compared with the manual method (which is not a gold standard method and can be different from one laboratory to another). However, the main advantage provided by the automatic peak fitting is the absence of intraoperator and interoperator variation and reproducibility. The second advantage of the automatic process is the normalisation of its preprocessing methods (MMA removal and baseline subtraction), removing the influence a human operator has in the spectrum preparation. Finally, the automatic method reduces drastically the time needed to analyse a spectrum (10 s is required for the automatic process, while the manual process takes 5 to 10 min per spectra).


## 3.1. Validation

The first step of the validation process is done by analysing the precision of the manual peak fitting method. Two types of precision are measured: (1) the intraoperator precision that measures the ability of one operator to obtain the same result multiple times on the same spectrum and (2) the interoperator precision that measures the ability to get the same result from two different operators on the same spectrum.The first operator (operator 1) is experienced with the bone FTIRM peak fitting, and is considered the reference. The second (operator 2) is an inexperienced operator who has been given the instructions to execute the peak fitting analysis.Five different bone spectra have been peak fitted with the manual method two times by each operator (cf. Table2). For the manual peak fitting, the root mean square of the coefficient of variation shows that results are operator-dependant.Table 2
Differences of the manual peak fitting analysis intraoperator, interoperator, maximum absolute error, and variation of results using the automatic method for each parameter, based on the analysis of 5 bone spectra. Operator 1 is an experienced operator, and operator 2 is an inexperienced operator to whom has been given the instructions to analyse the spectra.
Crystallinity Index
Mineral Maturity
Collagen Maturity
RMS of coefficient of variation (operator 1) (%)
0.7
1.4
0.7
RMS of coefficient of variation (operator 2) (%)
1.7
0.4
1.0
RMS of coefficient of variation (interoperator) (%)
2.0
1.2
0.7
Maximum absolute error
0.0013 cm
0.07
0.11
Variation of results using automatic method
No variation
No variation
No variation
RMS = root mean square.Then, we compared both automatic and manual methods based on the spectra extracted from the human samples and mouse sample. The Bland–Altman graph plots the means and the differences between automatic and manual methods measured by a human operator. For the human bones, Figure8 shows, at first sight, an almost null shift for every parameter. The ratio of valid results (i.e., results between the mean ± 1.96 ∗ SD) is 92% for the Crystallinity Index, 93% for the Mineral Maturity, and 92% for the Collagen Maturity.Figure 8
The Bland–Altman plot of the human bone spectra (manual method-automatic method values). The straight line is the mean of the differences. The dashed line is the mean of the differences plus or minus 1.96 ∗ SD, representing the limit of agreement. 92% of Crystallinity Index values, 93% of Mineral Maturity values, and 92% of Collagen Maturity values are within the dashed lines. Shifts of 0.0007 cm for the Crystallinity Index, 0.11 for the Mineral Maturity, and 0.17 for the Collagen Maturity are observed.
(a)
(b)
(c)For the mice tibia (Figure9), the shift is also almost null for every parameter. The ratio of valid results under the aforementioned hypothesis is 90% of valid results for the Crystallinity Index, 100% for the Mineral Maturity, and 100% for the Collagen Maturity.Figure 9
The Bland–Altman plot of the mice bone spectra (manual method-automatic method values). The straight line is the mean of the differences. The dashed line is the mean of the differences plus or minus 1.96 ∗ SD, representing the limit of agreement. 90% of Crystallinity Index values, 100% of Mineral Maturity values, and 100% of Collagen Maturity values are within the dashed lines. Shifts of −0.0009 cm for the Crystallinity Index, −0.01 for the Mineral Maturity, and 0.21 for the Collagen Maturity are observed.
(a)
(b)
(c)


## 3.2. Discussion

Gaussians have been long used to establish methods and peaks of interests in human bone, meaning changing the fitting function would end in a change of the peaks of interest. Thus, to maintain consistency with previous studies, Gaussians are kept for the fitting method. However, it should be noted that Lorentzians or Voigt functions would be more accurate regarding the quantum reality of the spectra [19].The development of the automatic method not only allowed a significant gain of time but also suppressed the interoperator and intraoperator variability.Currently, there is no gold standard to check the exactitude of the results obtained. After careful examination of the automatic peak fitting of these spectra, it has been found that the results obtained with the automatic method are valid compared with the analysis done by the expert using the manual method.The observed shifts are mainly due to the choice of boundaries within which the peak fitting method is applied. Indeed, there are no strict boundaries required for the manual peak fitting of theν4PO4 area. This induces the variation of the shift, as shown in Appendix B. This problem is solved by the automatic method that fixes those boundaries, chosen after careful analysis of several sets of boundaries.We could argue that it is not the most exact method in terms of localisation of the peaks, compared with the manual method (which is not a gold standard method and can be different from one laboratory to another). However, the main advantage provided by the automatic peak fitting is the absence of intraoperator and interoperator variation and reproducibility. The second advantage of the automatic process is the normalisation of its preprocessing methods (MMA removal and baseline subtraction), removing the influence a human operator has in the spectrum preparation. Finally, the automatic method reduces drastically the time needed to analyse a spectrum (10 s is required for the automatic process, while the manual process takes 5 to 10 min per spectra).


## 4. Conclusions

The algorithms described in this paper offer the possibility to analyse quickly an infrared spectrum of the thin bone section. This automatic method will allow the analysis of large sets of samples, which confer robust results.By using machine learning methods on large data of deconvolved FTIRM spectra, the automatic method could be improved. This could bring specificity of peak fitting parameters regarding the type of the sample analysed, thus improving the exactitude and reliability of the results.The software will be available under GNU GPL license.



---

*Source: 4131029-2018-12-23.xml*