QUANTIFICATION OF THE 3 D MORPHOLOGY OF THE BONE CELL NETWORK FROM SYNCHROTRON MICRO-CT IMAGES

In the context of bone diseases research, recent works have highlighted the crucial role of the osteocyte system. This system, hosted in the lacuno-canalicular network (LCN), plays a key role in the bone remodeling process. However, few data are available on the LCN due to the limitations of current microscopy techniques, and have mainly only been obtained from 2D histology sections. Here we present, for the first time, an automatic method to quantify the LCN in 3D from synchrotron radiation microtomography images. After segmentation of the LCN, two binary images are generated, one of lacunae (hosting the cell body) and one of canaliculi (small channels linking the lacunae). The binary image of lacunae is labeled, and for each object, lacunar descriptors are extracted after calculating the second order moments and the intrinsic volumes. Furthermore, we propose a specific method to quantify the ramification of canaliculi around each lacuna. To this aim, a signature of the numbers of canaliculi at different distances from the lacunar surface is estimated through the calculation of topological parameters. The proposed method was applied to the 3D SR micro-CT image of a human femoral mid-diaphysis bone sample. Statistical results are reported on 399 lacunae and their surrounding canaliculi.


INTRODUCTION
Bone diseases severely affect the quality of life around the world.According to a survey of International Osteoporosis Foundation, 1 of 3 women and 1 of 5 men with age over 50 suffer from osteoporosis, which is associated to bone fractures.Bone fragility remains only partially understood despite decades of research in this area.Recently, the crucial role of the osteocyte system at the cellular scale was highlighted (Bonewald, 2011).The osteocyte system, which is composed of osteocytes communicating through dendrites, is essential in bone mechanosensation and bone mechanotransduction and is supposed to play an important role in orchestrating bone adaptation (Burger and Klein-Nulend, 1999;Han et al., 2004;Gu et al., 2007).
The osteocyte system is hosted in the lacunocanalicular network (LCN).Although the LCN is raising increasing interest, it is a structure difficult to assess since it is deeply embedded in calcified bone matrix and it is composed of a huge number of nanometric structures.Lacunae are flat ellipsoidal voids housing the osteocyte bodies and with a density between 26000 to 90000/mm 3 (Cardoso et al., 2013).Canaliculi are narrow tunnels, enclosing the cell processes of the osteocytes with a reported diameter around 100-700 nm (You et al., 2004).The canaliculi connect lacunae and allow the circulation of interstitial fluid.Osteocytes are hypothesized to be stimulated by load-induced interstitial fluid displacement circulating in the LCN.Computational transport models have been proposed to estimate shear stresses (Weinbaum et al., 1994;Anderson et al., 2005;Anderson and Knothe Tate, 2008), but they are generally based on idealized lacuno-canalicular geometries.Therefore, the assessment of the LCN geometry is important to help understanding the strain-sensing mechanisms of the osteocyte network.
However, the 3D assessment of the LCN is still limited.Conventionally two-dimensional imaging techniques, such as light microscopy, atomic force microscopy (AFM), scanning electron microscopy (SEM) or transmission electron microscopy (TEM), have been proposed to visualize and quantify the LCN (Marotti et al., 1995;Kamioka et al., 2009;Lin and Xu, 2011;Sharma et al., 2012).Confocal microscopy can provide three-dimensional images, however, the limited penetration of light restricts the thickness of the visualized specimen and the spatial resolution is anisotropic (Sugawara et al., 2005).Recently, new 3D imaging techniques such as ptychography and serial FIB/SEM have been reported to image the LCN respectively at spatial resolutions of 65 nm (Dierolf et al., 2010) and 30 nm (Schneider et al., 2011).Due to the very high spatial resolution of these techniques, the field of view remains limited to a few cells, which restricts the possibility to obtain statistically significant results on the osteocyte network.
3D synchrotron radiation microtomography (SR micro-CT) has been postulated to be an ideal technique for visualizing the LCN (Müller, 2009).In previous works, we have demonstrated the feasibility of SR micro and nano-CT to image the LCN at the European Synchrotron Radiation Facility (ESRF) (Langer et al., 2012;Pacureanu et al., 2012).SR CT presents a number of advantages over conventional cone-beam CT.In particular, since the photons flux is several orders of magnitude higher than the one of conventional X-ray tubes, sub-micrometric spatial resolution can be reached while keeping good signal to noise ratio in relatively short time.With a 3D isotropic voxel size of 280 nm, this technique can provide a relatively large field of view (FOV) of about 600 3 µm 3 .Compared to the other 3D approaches, a few thousands of cells may be imaged in a single scan.
However, due to the novelty of these images and the complexity of the LCN, the current quantification tools cannot provide an automatic analysis of the LCN.Previous efforts have been made for the segmentation of LCN from the reconstructed images (Pacureanu et al., 2009(Pacureanu et al., , 2010(Pacureanu et al., , 2011;;Zuluaga et al., 2011).In order to quantify this network, techniques for extraction of dedicated parameters need to be designed.
Our aim is to propose new automated methods to quantify the LCN from such 3D SR micro-CT images.The next section briefly describes image acquisition and segmentation and presents the methods used to quantify the lacunae and canaliculi.In particular we present a new scheme to assess the distribution of canaliculi around each lacuna.The last section reports the validation of the method on a simple geometrical phantom, an isolated lacuna and finally results on a large volume encompassing 399 lacunae.

SAMPLE PREPARATION
The specimen was extracted from a woman (78 years old) femoral mid-diaphysis, obtained from a multi-organ collection.The sample was prepared following a series of treatments (Biobank, Presles en Brie, France), which underwent delipidation, elimination of medullary protein and sterilization.Then, the sample was stored at a room temperature until the imaging acquisition.The collection of sample followed the procedure of the Human Ethics Committee of the "Centre du don des Corps" at the University Rene Descartes (Paris, France) and is in accordance with legal clauses stated in the French Code of Public Health.

SYNCHROTRON RADIATION MICRO-CT IMAGING
Image acquisition was performed on the 3D parallel beam SR micro-CT setup developed at beamline ID19 at ESRF (Salomé et al., 1999).The pixel size was fixed to 280 nm providing a FOV of about (600 µm) 3 .The X-ray beam energy was set to 18 keV.To reduce the effect of radiation dose on the sample, it was necessary to optimize the detector setup (Pacureanu et al., 2012).A 4.5 µm thick LSO scintillator was used to transform the X-ray into the visible light which was detected with an E2V CCD camera.During the scan, the samples were rotated around a vertical axis with a total angle of view of 180° and 2000 projections were recorded.The total time for each scan was around 30 minutes.After scanning, the 3D image was reconstructed using a standard filtered backprojection algorithm.The reconstructed image has an isotropic voxel size of 280 nm.Fig. 1a illustrates a region of interest from a transversal slice of a reconstructed image (904 × 649 pixels).The organization of the lacuno-canalicular network concentrically distributed around the Haversian canal can be observed.

SEGMENTATION
The segmentation of the LCN from the micro-CT images is challenging due to the small size of canaliculi compared to the voxel size in the reconstructed image.Here, we followed a method introduced in a previous work to perform the segmentation (Pacureanu et al., 2010).First, a mask including the relevant parts of the bone tissue and excluding Haversian canals is computed.This was generated by using a large median filter (radius = 25) followed by a simple thresholding (see Fig. 1b).The segmentation of the LCN was performed through a variational region-growing method (Pacureanu et al. 2010;Revol-Muller et al., 2012) based on an energy functional combining grey level information from the original image and shape information extracted via a 3D tube enhancement filter based on the eigenvalues of the Hessian (Sato et al., 1998;Pacureanu et al., 2013).Subsequently, the lacunae and canaliculi are separated by a 3D opening operation, with a 3 × 3 × 3 ball-shaped structuring element.After this step, two binarized images, of the lacunae and the canaliculi are available.Fig. 1c illustrates a composite image of the segmented LCN.
We can formalize the segmented LCN as the union of each lacuna L n with its set of canaliculi C n : ( ) where M is the total number of lacunae in the image.

QUANTIFICATION METHOD Labeling
In order to quantify each lacuna L n with its set of canaliculi C n , the segmented lacunae are labeled by using a connected component analysis.Here, a fast labeling method proposed by Hoshen and Kopelmanm (1976) was implemented in 3D.This method allows labeling of all objects by only scanning the image twice in a raster pattern, i.e. from top to bottom and bottom to top.

Quantification of lacunae
After labeling, the total number of the lacunae (N.Lc) and the lacunar density with respect to the bone volume (N.Lc/BV) and the tissue volume (N.Lc/TV) can be easily calculated.The bone volume (BV), the tissue volume (TV) and the canal volume fraction (Ca.V/TV) were quantified from the mask image.Then, a number of descriptors were extracted from each labeled lacuna L n as described below.

Moment based descriptors
Considering the ellipsoidal shape of the osteocyte lacunae, we used a 3D moment-based approach to fit each lacuna to an ellipsoid (Denis et al., 2008;Flusser et al., 2009).Let μ pqr (L n ) be the second-order central moments of L n , with p+q+r=2 for and M(L n ) be the corresponding 3 × 3 second order central moment matrix.
We denote λ k (k = 1,2,3) the eigenvalues of M(L n ) sorted by decreasing value.The half axes of the fitting ellipsoid having the same moment matrix are given by: ( ) where Lc.V(L n ) is the lacunar volume, which equals to the total mass of the object (zero-th order moment).
From the fitting ellipsoid, the osteocyte lacunar length (Lc.L 1 (L n ) = 2a 1 ), width (Lc.L 2 (L n ) = 2a 2 ), and depth (Lc.L 3 (L n ) = 2a 3 ) were calculated.In addition, the lacunar anisotropy was quantified as The ratio between the volume of the fitting ellipsoid and the actual one is given by: Besides, we also calculated the lacunar volume ratio with respect to the bone volume (Lc.V/BV) and the tissue volume (Lc.V/TV).Finally, the distance between nearest lacunae (Lc-Lc) was estimated by computing for each lacunae the minimum distance between its gravity center to that of the closest osteocyte lacuna.

Intrinsic volumes
In addition, we calculated the intrinsic volumes to further quantify the surface area Lc.S(L n ), the Euler number Lc.χ(L n ) and the structure model index (SMI) Lc.SMI(L n ) of the lacuna L n (Ohser and Schladitz, 2009).The intrinsic volumes are invariant geometric functions serving as a basis of object features.In our three-dimensional image with the orthorhombic cubic lattice of spacing Δ, the four intrinsic volumes , respectively representing the Euler number, integral of mean curvature, surface area and volume of the lacuna, can be efficiently calculated from the 3 , expressed as: where l is an index spanning local configurations of voxels, is the number of configuration l in the set L n and ) (k l v is the weight coefficient of the local configuration.The surface area Lc.S(L n ), the integral of the mean curvature Lc.M(L n ), and the Euler number Lc.χ(L n ) of the lacuna L n can be estimated as: , and the weights ) 3 ( l v for the Euler number for each configuration are given in (Ohser et al., 2009).
In addition, we also calculated the structure model index of each lacuna (Lc.SMI(L n )) from the intrinsic volumes (Ohser and Schladitz, 2009), expressed as: It characterizes the lacunar shape with values of 0, for a pure plate, 3 for a rod and 4 for a sphere (Hildebrand and Rüegsegger, 1997).

Quantification of ramification of canaliculi
We propose a new method to characterize the 3D ramification structure of canaliculi issued from a lacu-na.Considering that the canaliculi are branching after radiating from the surface of lacunae, our goal is to give a signature of the number of canaliculi per lacuna at different distances from the lacuna boundary.
This method is based on the use of the 3D Euler number.Let us recall that the Euler number of an object in a three-dimensional space can be expressed as a function of the Betti numbers β 0 , β 1 , and β 2 as: where β 0 , β 1 , β 2 represent respectively the number of connected components, the number of tunnels and the number of cavities in the object (Odgaard, 1997).
In practice, when the object is represented as a discrete set of voxels, the 3D Euler number can be calculated as (Serra, 1982;Toriwaki et al., 2002): where n 0 , n 1 , n 2 , and n 3 are respectively the number of vertices, edges, faces, and voxels.
To calculate the number of canaliculi per lacuna at different distances from the lacunar surface, first we extract the bounding surface of the lacuna using morphologic derivative operations.The bounding surface ∂ r L n can be formulated as the subtraction between the r + 1 and r times dilated lacuna: where ⊕ denotes dilation.This bounding surface will be used to count the number of canaliculi at the distance r, from the surface of the lacuna.Each bounding surface ∂ r L n has a single connected component with one cavity inside.
Second, a special bounding surface with holes is generated as: After that, the Euler number of the bounding surface with holes Finally, the number of canaliculi per lacuna ) )( ( .r L NCa Lc n at a given distance r can be calculated as: The dilation parameter reflects the distance at which the count is performed.Thus for each lacuna, gives a signature of the distribution of canaliculi around the lacuna.The r max was chosen smaller than the mean inter lacunar dis-tance Lc-Lc.The evolution of the number of cana-liculi with r reflects the ramification of the canaliculi. To further characterize the branching of canaliculi, we also calculated the ratio between the maximum and minimum number of canaliculi per lacuna at different distances from the surface of lacunae.Besides, we estimated an average length of the primary canaliculi for a given lacuna r p (L n ).The calculation was done under the assumption that the probability of any canaliculi to bifurcate is uniform in a certain range of distance.When considering only the first bifurcation, the number of canaliculi with respect to the distance (r) to the surface of lacunae increases linearly.Therefore, the relation between the number of canaliculi versus the distance was fit by a straight line Lc.NCa(L n )(r) ~ a(L n )r + b(L n ).Then, r p (L n ) is estimated as the length for which 1.5 times the initial number of canaliculi Lc.NCa(L n )(r 0 ) and by correcting the distance step ∆r/2: Finally, the descriptive statistics of each parameter were calculated on a whole population of lacunae in a given bone sample.

VALIDATION OF THE CANALICULI COUNTING METHOD
The quantification of numbers of canaliculi was validated on two test volumes.The first one is a simple geometrical phantom constituted of three ellipses each with six cylinders, simulating three lacunae interconnected by canaliculi (cf.Fig. 3a).The volume size is 64 3 voxels.The second is the segmentation of one isolated lacuna with all canaliculi interconnected, shown in Fig. 3b.The volume size is 149 × 149 × 85 voxels.
The application of the method to the first simple phantom provided 6 canaliculi for each ellipse as expected by construction.Concerning the second experimental image of an isolated lacuna, the validation was carried out by manual check and by a test program.The test program consisted in applying a connected component analysis after removing the lacuna.The method was applied with two choices of the dilation parameter r.For r equals 1 (respectively 15), the number of canaliculi found was 22 (respectively 32).The comparison with the test program in these two cases (Fig. 3c,d) shows that the correct number of canaliculi was obtained.The difference between the two results reveals that canaliculi are branched as can be observed in the 3D renderings.Thus, the dilation parameter is useful to put in evidence the ramification of canaliculi.

APPLICATION TO SR MICRO-CT IMAGE
The proposed method was applied to a large SRmicro-CT image acquired at the ESRF at a spatial resolution of 280 nm.A Volume of Interest (VOI) made of 904×649×998 voxels was extracted from a (2048) 3 reconstructed volume.The VOI includes an osteon with hundreds of lacunae Fig. 1a.In the segmented osteocyte lacunae image, unacceptable objects such as micro-cracks, ring artifacts, and some irregu-larly shaped fragments can be extracted as well.From the literature, it is known that the size and shape of lacunae resides in a certain range (Mc-Creadie et al., 2004;Dong et al., 2014).Thus, the segmented lacunae were filtered by thresholding some descriptors as follows: 50 µm 3 < Lc.V < 1000 µm 3 , Lc.L1/Lc.L2 < 5, and 0 ≤ Lc.χ ≤ 2. In addition, we also removed the outliers, by excluding the lacunae for which the canalicular regression slope was smaller than 0.5 or higher than 50.The segmented image included 399 lacunae and their associated canaliculi.Fig. 4 illustrates the labeled lacunae in different colors.We report the bone tissue histomorphometric parameters in Table 1.Besides, the statistical results on the lacunar descriptors with the mean, standard deviation, minimum and maximum are reported in Table 2.The average distance between nearest lacunae was calculated and was found to be 23.2 µm.
The mean and standard deviation of lacunar volume and surface was 216.4 ± 84.7 µm 3 and 238.9 ± 66.4 µm 3 .The average lengths of the three axes of lacuna are 15.2 µm, 7.8 µm and 4.0 µm, and length ratio between the three axes is about 4:2:1.The mean of lacunar SMI is 3.2.
To quantify the numbers of canaliculi, we used 8 different dilation parameters (r), from 5 to 40 with step of 5 corresponding to distances from 1.4 µm to 11.2 µm with step of 1.4 µm.Boundary lacunae were excluded from analysis, due to incomplete LCN.The statistical results are reported in Fig. 5.
The average number of canaliculi per lacuna was found to increase between 41.5 and 139.1 with the distance r.The linear regression explains 98.0% of the variability of the number of canaliculi around their mean (see Fig. 5).The average length of the primary canaliculi is about 5.6 µm.The median of ratio between the maximum and minimum number of canaliculi per lacuna is 3.4.Fig. 6 illustrates 3 different lacunae with different number of canaliculi.Fig. 6a shows a lacuna with large number of 98 canaliculi, Fig. 6b, the one with most frequent number of 65 canaliculi, and Fig. 6c, a lacuna with a small number of 49 canaliculi.The three ramification patterns calculated for each of these lacunae are also shown in Fig. 6d.

DISCUSSION
In this paper, we proposed a technique to automatically calculate LCN descriptors from a segmented 3D SR micro-CT image.Both lacunae and canaliculi were quantified.Considering the lacunae, we calculated both conventional shape parameters and parameters based on ellipsoidal fit.An original method to assess the ramification pattern of the canaliculi was proposed based on mathematical morphology operators and the Euler number.It consists in counting the numbers of canaliculi per lacunae at different distances.Conversely to a method that would try to identify each node of the canaliculi network, this method has the advantage to be more flexible to imperfections in segmentation which are unavoidable with the limited spatial resolution compared to the canaliculi size.
The estimation of the numbers of canaliculi per lacunae was checked on a simple phantom and a complete lacuna with fully interconnected canaliculi.Then, this method was applied to a larger experimental image of human bone tissue.This is the first time that the numbers of canaliculi at different distances are measured in 3D.We have extracted these measurements on several hundreds of cells in human femoral bone.Traditional investigations are based on 2D histological sections analyzed manually.Compared to the extrapolation of 2D measurements to 3D, the present method does not require any assumption on idealized lacuno-canalicular models.Thus this method is expected to yield unbiased parameters.Since the method is fully automated, it can be applied to bone samples within large field of view including several osteons.
The canaliculi were counted in 240 lacunae, after excluding boundary lacunae.Keeping all lacunae would provide biased results due to truncation of LCN.The average number of canaliculi near the surface of lacunae is in agreement with the results published by Beno (Beno et al., 2006) who reported an estimated 3D number of canaliculi of 41 in human bone.The low and high estimates are also consistent with Beno's work.In a more recent work, the number of canaliculi per lacuna was calculated from confocal microscopy images on a limited depth of 25 µm (Sharma et al., 2012).The authors reported about 85 primary and 387 secondary canaliculi per lacuna on rats' tibial metaphysis.These higher numbers can be explained by differences in method, species and bone site.In our method, the various dilation parameters permitted to put in evidence the increased branching of canaliculi with the radial distance to the lacuna.Our results clearly indicate an increase of the number of canaliculi with the distance to the lacunae.In addition, for most lacunae, the increase could be explained by a linear regression between the number of canaliculi and the distance suggesting the bifurcation of the canaliculi.Table 1.Histomorphometric parameters of the bone tissue from the human femoral sample.
Lc.TV BV TV BV/TV Ca.V/TV Lc.TV/BV Lc.TV/TV N.Lc/BV N.Lc/TV N.Lc (mm 3 ) (mm 3 ) (mm 3 ) (%) (%) (%) (%) (mm -3 ) (mm -3 ) 399 8.63 × 10 -5 12.39 × 10 -3 12.85 × 10 -3 96.4% 3.6% 0.70% 0.67% 32200 31042 N.Lc -number of osteocyte lacunae; Lc.TV -total lacunae volume (mm 3 ); BV -bone volume (mm 3 ); TV -tissue volume (mm 3 ); BV/TV -bone volume fraction (%); Ca.V/TV -canal volume fraction or bone porosity (%); Lc.TV/BV and Lc.TV/TV-lacunar volume density (mm -3 ); N.Lc/BV and N.Lc/TV-lacunar number density (mm -3 ) The results strongly rely on image segmentation which itself is dependent on the quality of the acquired images.Further advances in optimizing the imaging technique may improve image quality, while progress in the segmentation methods would reduce noise and improve the connectivity of the lacunocanalicular network.Although the numbers of canaliculi that we found are realistic, it is difficult to evaluate the method due to the lack of comparison data and due to the complexity of the network.These preliminary observations need to be confirmed by further works.
The method has the advantages of producing 3D measurements in an automated and model independent manner.We believe that this can contribute to improving our knowledge on the LCN.Furthermore, future work will also be pursued to extract additional characteristics of the LCN, such as the canaliculi length and the canaliculi density per each lacuna.
Fig. 1.(a) A slice of reconstructed image of SR micro-CT (voxel size 280 nm, FOV: 904 × 649 pixels), including a Haversian canal, lacunae and canaliculi.(b) A slice of mask image which excludes the region of Haversian canal.(bone tissue in white and porosity in black) (c) The corresponding segmented lacunae (in yellow) and canaliculi (in green).

Fig. 2 .
Fig. 2. (a) Segmented image with 3 flat elliptic-like lacunae and many tube-like canaliculi.(b) Labeled lacunae after the connected component analysis.Different colors show the different labels.The 3D renderings were generated with the software Avizo®.
Fig. 4. 3D rendering of the labeled lacunae (displayed in different colors) concentric distributed around the osteon.(a) top view (b) side view.

Fig. 5 .
Fig. 5.The mean and standard deviation of the number of canaliculi counted at the eight different distances from the surface of lacunae (n = 240; boundary lacunae were excluded from analysis).

Fig. 6 .
Fig. 6. 3D renderings of three isolated lacunae cropped from the segmented image.(a) A lacuna with 98 canaliculi counted at 5.6 µm to surface of lacuna.(b) A lacuna with 65 canaliculi counted at 5.6 µm to surface of lacuna.(c) A lacuna on the boundary of the osteon with 49 canaliculi counted at 5.6 µm to surface of lacuna.(d) Plots of the number of canaliculi calculated at eight different distances from the surface of the three selected lacunae.

Table 2 .
Osteocyte lacunar descriptors from the human femoral sample.