Brillouin scattering from biomedical samples: the challenge of heterogeneity

Brillouin light scattering (BLS), a non-destructive and non-contact technique, offers a powerful tool for probing the micromechanical properties of biological tissues. However, the inherent heterogeneity of biological tissues can pose significant challenges in interpreting BLS spectra. In this study, we introduce a novel method that harnesses the intensity information within a single BLS spectrum to directly estimate the Voigt average of the longitudinal modulus. Additionally, we use a method to determine the ratio of the squared Pockels coefficients for photoelastically heterogeneous samples, based on global analysis of a 2D BLS map. This method is shown to effectively determine the photoelastic ratio of soft and hard components of human bone tissues, enabling the calculation of the average elastic moduli. Furthermore, it has the remarkable ability to generate maps of the filling factor of the scattering volume, shedding valuable light on the intricate structure and topography of rough surfaces under BLS mapping.


Introduction
Comprehending the mechanical properties of biological materials is of paramount importance for advancements in healthcare, biotechnology, and numerous other fields.However, this task presents significant challenges due to the intricate and dynamic nature of these materials.To adequately capture their mechanical properties across various scales, a combination of techniques is often required.In recent years, Brillouin light scattering (BLS) has emerged as a powerful tool for viscoelastic characterization at the microscale, offering a non-destructive and non-contact approach to probe the mechanical behavior of biological materials at a microscopic level [1][2][3][4][5][6][7].Advances in BLS technology have led to significant improvements in the speed of data acquisition, enabling rapid Brillouin spectral imaging [8][9][10][11][12][13].This breakthrough has fuelled the growing demand for BLS hyperspectral viscoelastic imaging of biomedical samples.In fact, BLS imaging can accelerate research in areas such as biomechanics [14,15], mechanobiology [16,17], and tissue engineering [18,19] by providing access to valuable information about the biomechanical properties of cells and tissues [20][21][22].It is evident that enhancing the reliability of spectral interpretation can significantly expand the clinical utility of BLS imaging, paving the way for innovative diagnostic tools and for evaluating tissue health and disease [22][23][24][25].
However, the interpretation of Brillouin spectra remains a challenge, particularly for complex biological tissues with intricate microstructures.Unlike homogeneous condensed matter, biological tissues exhibit a high degree of heterogeneity, often displaying hierarchical structures and viscoelastic behavior, characterized by complex elastic moduli that vary depending on position and on timescale of the mechanical perturbation.Bone tissue stands as an exemplar of hierarchical spatial structures, exhibiting both sub-and super-micrometric heterogeneities that give rise to multi-peak spectral features [26][27][28][29].These intricate patterns hold valuable information about the nature and shape of the alternating soft and hard structures that define the tissue's architecture.However, the sheer volume of information packed into a single spectrum presents novel challenges in interpreting the spectral features, unlike those encountered in traditional Brillouin scattering studies of homogeneous condensed matter.We are now confronted with the complexities of, e.g.heterogeneous broadening of Brillouin peaks due to longitudinal acoustic modes propagation within different structures, possibly exceeding the region lit by the incoming laser beam [30] and with the need of estimating average elastic moduli in presence of both elastic and photo-elastic heterogeneities.Here, we give some new hints and summarize recent advances in these key topics, a contribution to unlocking the full potential of micro-Brillouin mapping in biological matter.Our experimental setup is described in [31].

Heterogeneity in space and heterogeneous broadening
Biomedical samples often exhibit spatial heterogeneity at the micrometer scale, resulting in a single spectrum containing Brillouin peaks generated from distinct subregions.
A primary concern is that micrometer-scale heterogeneities can lead to multiple scattering of the probing light.This phenomenon is particularly pronounced in deeper regions of the sample, where it contributes to a progressive increase in spurious signals arising from randomly distributed scattering wavevectors.A new experimental strategy, the so-called Polarization Gated Brillouin Spectroscopy PG-BS, was recently proposed to address this issue [32].In the following, we focus on the simpler case where scattering occurs under negligible multiple scattering conditions, such as in scattering volumes near the sample surface.
Whenever possible, distinguishing the contributions from each component is advantageous when analysing Brillouin peaks generated by different subregions of the sample.This approach provides valuable insights into the viscoelastic properties of the individual components [33,34], but it is not always possible.For example, when the dimension of the heterogeneities are smaller or comparable with the acoustic wavelength and the acoustic mismatch between the heterogeneities and the surrounding media is low, their presence is not directly detectable as a different peak in the spectrum.A single acoustic mode can actually span across the different structures directly averaging the elastic properties of the medium.This phenomenon has been studied in heterogeneous amorphous systems, where heterogeneities have been shown to produce an additional attenuation process of the acoustic waves-the so-called static contribution-which adds to the dynamic ones [35,36].When the dimensions of the heterogeneities become comparable to the acoustic wavelength and the acoustic impedance mismatch between the heterogeneities and the surrounding medium is sufficiently high, localized acoustic modes can be confined within these structures, resulting in a multipeaked spectrum.The frequency positions of these peaks provide information about the elastic properties of the heterogeneities, but this information can be fully extracted only in conjunction with knowledge of their characteristic size [37].The emergence of distinct spectral peaks, whose positions are directly linked to the sound velocities of the individual structures, requires the dimensions of these heterogeneities to exceed the coherence length scale of acoustic modes.This length scale represents the characteristic size over which density fluctuations, the primary quantities probed by BLS, exhibit correlation.An example is reported in figure 1.
When this condition is met, a multi-peak analysis can be employed to derive the mechanical properties of the constituent structures.However, even in this case, deconvolution is not always possible, particularly when the spectrum exhibits multiple closely spaced peaks that are challenging to differentiate.This situation arises in cases of 'heterogeneous broadening,' where the spectral bandwidth cannot be attributed solely to the viscous properties of the sample.In such instances, providing an average value for the elastic modulus is a reasonable outcome of the analysis.The average elastic modulus is usually determined by fitting the spectrum with a single broad peak or by calculating the first spectral moment.Both approaches yield an average frequency ⟨ω 0 ⟩, whose squared value is then used to estimate the average elastic modulus, assuming ⟨M⟩ = ρ/q 2 ⟨ω 0 ⟩ 2 , where ρ is the mass density and q is the wavevector of the longitudinal acoustic modes that, in back scattering experiments, is given by q = 2nk i with n the refractive index of the medium and k i the wavevector of laser light in vacuum.These methods are suitable for relatively narrow distributions of Brillouin peaks but become increasingly inaccurate as the distribution broadens.The average elastic modulus should ideally be modeled as a weighted average of the moduli of the different components, based on the rule of mixtures.Here we propose the Voigt average, which has been shown to better fit hydrated collagen elastic moduli than other elastic models [39].For an N-component system, the Voight average is defined as is the elastic modulus of one component with Brillouin peak at ω 0i and v i is its volume fraction.In case where the ratio ρ/q 2 is not strongly dependent on the component, as observed in some biological systems [40,41], we can write: It can be seen that the average elastic modulus is related to the average of the squared frequencies ⟨ ω 2 0 ⟩ rather than the square of the averaged frequencies ⟨ω 0 ⟩ 2 .To address this issue, we propose a rapid and effective method to estimate the Voigt average of the longitudinal modulus directly from the Brillouin spectrum of an heterogeneous sample, utilizing the intensity of the spectrum as a proxy for the volume fractions of the sample components.Under specific conditions, the weighted average of ω 2 0i can be calculated and used to estimate the average elastic modulus of the sample.These conditions are: (i) The Voigt method provides a reliable estimate of the average elastic modulus; (ii) The major contribution to the variation of ω 2 0i arises from the elastic modulus rather than the ratio ρ/q 2 ; (iii) The photoelastic constants (Pockels coefficients) of the averaged subcomponents of the sample are not significantly different from one another.If the aforementioned conditions are satisfied, the weighted average of ω 2 0i can be determined using the intensity of the spectrum as a weighting factor, following the subtraction of any background contribution: where ω j and I j are the frequency and the intensity of a single point of the spectrum.The value of ⟨M⟩ can be estimated as: It is worth noting that this method can be employed under both conditions of heterogeneous broadening of a Brillouin peak and in cases where separate peaks are simultaneously present in the spectrum.In fact, this applicability extends to all situations where the average elastic modulus is of interest, rather than the individual elastic moduli of the components.
To estimate the order of magnitude of the improvement obtained with calculating we calculate the averages of two DHO functions with same intensity separated by a frequency difference of ∆ω = ω 02 − ω 01 (figure 2).For our purpose, the two peaks can be considered either as the higher (ω 02 ) and lower (ω 01 ) frequency limits of a given distribution, or as two separate components of the sample to be averaged.Numerical integration of two DHO functions was performed within the frequency range of 4-10 GHz.The DHO functions were selected with linewidths Γ /2π = 0.1 GHz and average frequency ⟨ω⟩/2π = 6.5 GHz.Note that the relative difference between the expected value of ω 2 0 and the value calculated using equation ( 2) is found to be less than 0.1% for ∆ω/⟨ω⟩ up to 0.5 (filled squares in the inset of figure 2).This method markedly enhances the reliability of the estimated value of the elastic modulus compared to existing strategies that rely on calculating ⟨ω 0 ⟩ 2 .Indeed, a simple calculation shows that the relative error in substituting ⟨ ω 2 0 ⟩ with ⟨ω 0 ⟩ 2 is given by: This behavior is illustrated by the solid line in the inset of figure 2. The same graph also shows the values of ( ⟨ /⟨ω 0 ⟩ 2 obtained through the numerical calculation of ⟨ω 0 ⟩ from the first spectral moment of the spectrum (circles).The residual error values for ∆ω/ ⟨ω 0 ⟩ approaching zero can be attributed to the asymmetric nature of the DHO function.
It is important to note that the rapid and efficient averaging method of equations ( 2) and ( 3) is only suitable for Brillouin peaks that arise from longitudinal acoustic modes in subsamples with comparable photoelastic constants.For instance, when considering Brillouin scattering from bones, this approach may be valid within the 'soft' or 'hard' regions of the spectrum (see figure 1).However, averaging the entire spectrum would not be correct due to the significantly different photoelastic constants of the soft and hard bone components.To address this limitation, we recently introduced a method for estimating the ratio of the squared Pockels coefficients by analyzing the entire set of spectra gathered for a map [28], as described in the following section.

Photoelastic heterogeneity and filling factor
In a biphasic system composed of coexisting soft and hard components, as illustrated in figure 1, the average longitudinal modulus at the measured point can be obtained from equation (3) and the average ω 2 0 from equation: where are the average values of the hard and soft components obtained by equation ( 2).The volume fraction of the hard component v H can be determined using the following relationship: Where I H and I S are the integrated intensities of the Brillouin peaks and r is the ratio of the squared Pockel coefficient for the hard and soft portions of the sample.Its value can be obtained from the ratio between the maximum measurable intensities r = I HM /I SM , where I HM and I SM , ideally correspond to the intensity of the single components filling the whole scattering volume.
The Brillouin spectra used to test our methodology were from a small portion of a human femoral diaphysis and were already published in [38].The sample was obtained from a cadaver donor by the Musculoskeletal Tissue Bank of the IRCCS Istituto Ortopedico Rizzoli (Bologna, Italy; EU TE code: IT000096).According to Italian legislation and to the guidelines of the National Transplant Centre, as the sample was not suitable for transplantation and considered as waste material, it was only used for validation purposes.Further details on the structure of the diaphysis (cortical tissue), sample pre-treatment and fixation are provided in [26].Examples of spectra obtained at room temperature from the raster scan of a 185 × 96 µm section are shown in the insets of figure 3.Each spectrum of the map was acquired for 50 s.To  It is interesting to note that the fluctuations in the total integrated intensity of the spectrum evidenced by the scattered distribution of I H vs. I S values in figure 3 can give precious information on the filling factor of the scattering volume, which gives an estimate of the topology of the surface.In fact, the filling factor can be calculated as The map of the calculated values of R, reported in figure 4(e), evidences a rough surface, likely due to both the mineralized nature of bone tissue and the surgical resection.The Haversian canal, which contains the blood vessel in the center of the map, and the lacunae, i.e. the niches where osteocytes are located, are well recognizable due to the lower values of the filling factor.

Conclusions
In biological tissues, the presence of multiple components with varying elastic properties can make it difficult to accurately determine the average elastic modulus using traditional BLS methods.In this work, we have introduced a novel technique that utilizes the intensity of the BLS spectrum to estimate the Voigt average of the longitudinal modulus directly from the sample's Brillouin spectrum.This method, based on the intensity as a proxy for the volume fractions of the sample components, outperforms existing methods that rely on averaging the squared frequencies of the Brillouin peaks.Additionally, we have applied a recently developed approach to estimate the ratio of the squared Pockels coefficients for photoelastically heterogeneous samples.This method, which analyses the entire set of spectra gathered for a 2D BLS map, has been shown to be effective in determining the photoelastic ratio of soft and hard components of human bone tissues.Moreover, we have demonstrated its capability to generate maps of the scattering volume filling factor, providing valuable insights into the intricacies of the surface topography as observed through BLS mapping.We anticipate that this work will foster the advancement of BioBrillouin as an increasingly rapid, robust, and dependable technique.

Figure 1 .
Figure 1.(a) Schematic of the setup for Brillouin light scattering measurements.A high resolution Tandem Fabry-Perot interferometer is employed to analyze the spectral features of light back-scattered from bone tissue.(b) Graphical representation of the spatial heterogeneity inherent to bone tissue, arising from its micrometric constituents.Acoustic moded in each region can scatter light at different frequencies; (c) Brillouin scattering spectra acquired from a section of human femoral diaphysis[38], revealing two distinct components-P SOFT and PHARD-mainly attributed to the osteocytes (bone cells) and the mineralized extracellular matrix, respectively.Reproduced from[1].CC BY 4.0.Reproduced from[26].CC BY 4.0.

Figure 3 .
Figure 3. Graph depicting the integrated intensity of the hard component (IH) as a function of the soft component (IS) for each spectrum recorded on the map.Insets show in blue, green, and magenta, characteristic spectra relative to the corresponding points in the graph.The slope of the red line gives the value of the photoelastic ratio r to be used in equation (6).The dashed blue line and the values of I * H , I f Hand IHM determine the filling factor R at the point in the map corresponding to the green spectrum, using equation(7).Both in individual points of a spectrum and in integrated intensities, the error bars are given by the square root of the photon count[42].

Figure 4 . 2 0ω 2 0 ⟩
Figure 4. 2D maps of a 185 × 96 µm section of the human femoral diaphysis here analyzed, showing the relative intensity of the Brillouin peaks (a) and (c) and the corresponding means of the squares of the Brillouin frequencies (b) and (d) (here ν = ω/2π) of the soft and hard components, respectively; (f) the mean of soft and hard components obtained through equation (5); (e) the filling factor R from equation (7); (g) a graphical representation of the analyzed with its main components.