Estimation of Photon Path Length and Penetration Depth in Articular Cartilage Zonal Architecture Over the Therapeutic Window

Objective: The key characteristics of light propagation are the average penetration depth, average maximum penetration depth, average maximum lateral spread, and average path length of photons. These parameters depend on tissue optical properties and, thus, on the pathological state of the tissue. Hence, they could provide diagnostic information on tissue integrity. This study investigates these parameters for articular cartilage which has a complex structure. Methods: We utilize Monte Carlo simulation to simulate photon trajectories in articular cartilage and estimate the average values of the light propagation parameters (penetration depth, maximum penetration depth, maximum lateral spread, and path length) in the spectral band of 400–1400 nm based on the optical properties of articular cartilage zonal layers and bulk tissue. Results: Our findings suggest that photons in the visible band probe a localized small volume of articular cartilage superficial and middle zones, while those in the NIR band penetrate deeper into the tissue and have larger lateral spread. In addition, we demonstrate that a simple model of articular cartilage tissue, based on the optical properties of the bulk tissue, is capable to provide an accurate description of the light-tissue interaction in articular cartilage. Conclusion: The results indicate that as the photons in the spectral band of 400–1400 nm can reach the full depth of articular cartilage matrix, they can provide viable information on its pathological state. Therefore, diffuse optical spectroscopy holds significant importance for objectively assessing articular cartilage health. Significance: In this study, for the first time, we estimate the light propagation parameters in articular cartilage.


I. INTRODUCTION
L IGHT propagation in biological tissues is a complex pro- cess that holds significant importance in various fields of biophotonics, including diagnostics and therapeutics [1], [2].When light interacts with a tissue, its behavior is governed by intrinsic optical properties, such as absorption coefficient (μ a ), scattering coefficient (μ s ), scattering anisotropy factor (g), and refractive index (n) [2].These parameters determine whether photons are absorbed within the tissue or scattered in new directions, either towards the deeper regions or reflected at the tissue surface [1].Therefore, the characterization of light propagation within biological tissues relies on their optical properties, providing valuable information about their pathological status [3].As tissue optical properties can vary due to microstructure or pathological changes, understanding the light propagation characteristics of biological tissues can reveal insights into their structural and biochemical properties.The key characteristics of light propagation are the averaged penetration depth (P D avg ), averaged maximum penetration depth (P D max ), averaged maximum lateral spread (LS max ) and averaged path length (P L avg ) of photons.P D avg is the point in tissue depth where the weight of propagating photons falls to 63.21% of their incident weight, P D max represents the averaged maximum depth photons can reach, LS max is the averaged maximum reach of the photons on the transverse plane, and P L avg is the averaged distance photons travel inside the tissue volume [4], [5].
An example of biological tissue with complex structural and chemical properties is articular cartilage.It comprises an extracellular matrix (ECM) and chondrocytes -the cellular component of the tissue.The ECM mainly comprises collagen fibers, proteoglycan macromolecules, and water [6].Articular cartilage exhibits a layered architecture with distinct zones: superficial zone (SZ), middle zone (MZ), and deep zone (DZ) [7].The collagen fibers in the SZ align parallel to the articular surface, followed by a skewed alignment in the MZ and ultimately a perpendicular alignment to the articular surface in the DZ [8].Given the unique architecture of articular cartilage ECM, we aim to investigate how light propagation in bulk articular cartilage tissue is influenced by its zonal architecture.
It is known that articular cartilage zonal architecture influences its optical properties [9], [10].This study investigates the parametric characterization of light propagation in articular cartilage.Since the optical properties of articular cartilage and its zonal layers show wavelength-dependency, we hypothesize that P D avg , P D max , LS avg , and P L avg of photons across the visible and near-infrared (NIR) spectral regions would differ.As a result, we expect nonuniform light propagation properties across these spectral bands.To test this hypothesis, we used Monte Carlo (MC) method [1], [2], [11], [12], [13] to simulate photon trajectories and estimate their P D avg , P D max , LS avg , and P L avg in the spectral band of 400-1400 nm.By comparing these values with the thickness of the articular cartilage tissue, we aim to provide insight on how visible and NIR light can probe the zonal layers of articular cartilage ECM.
This study addresses the gap in knowledge regarding the characteristics of light propagation in articular cartilage, with potential implications for diagnostics and therapeutics.The findings provide wavelength-dependent information about P D avg , P D max , LS avg , and P L avg of photons within the spectral regions used for diagnostic and therapeutic purposes.The precise knowledge of these parameters is crucial for the advancement of medical device instrumentation, especially in the field of orthopedic diagnostics.The joint surgeries require accurate and real-time intraoperative information to expedite objective decision making for treatment planning.Non-destructive optical technologies such as near-infrared spectroscopy have shown great potential to bridge this gap.On the other hand, P D avg , P D max , LS max , and P L avg play a pivotal role in the design and optimization of optical-based medical devices tailored for orthopedic procedures.They provide crucial information about how much the diagnostic light penetrates the tissue axially and laterally, detailing the subsurface volume of the tissue that has been probed by light.Therefore, these parameters are essential for the design and instrumentation of diagnostic devices for characterizing the subsurface volume of musculoskeletal tissues during orthopedic surgeries, aiding surgeons in making informed decisions during surgeries.Thus, understanding the behavior of light in biological tissues is not only essential for the advancement of medical technologies, such as photoacoustic imaging and diffuse optical tomography but also holds promise for enhancing the accuracy and effectiveness of intraoperative joint procedures [14], [15], [16].

II. METHODS
In order to rigorously test our hypothesis, we employed a multifaceted methodology designed to comprehensively assess the optical and structural characteristics of articular cartilage.Firstly, we harvested two sets of articular cartilage samples sourced from distinct anatomical locations within bovine knees.
One set underwent microCT imaging to obtain accurate physical dimensions, essential for subsequent MC simulations.Subsequently we removed the end bones of the samples and used the remaining bulk articular cartilage part of the samples, in conjunction with an integrating sphere setup, to measure their optical properties, specifically μ a and μ s .Following optical properties measurements, the samples underwent histological processing to obtain thin axial histology sections of articular cartilage, crucial for polarized light microscopy (PLM) imaging.We employed PLM imaging to quantify angular orientation of collagen fibers in the tissue in order to objectively classify the tissue matrix into three layers of SZ, MZ, and DZ and understand how much of the tissue volume is occupied by each layer.
Simultaneously, the second set of harvested samples was utilized to estimate the μ s of individual articular cartilage zones.We first subjected the samples to cryosectioning, involving the acquisition of thin lateral sections from articular cartilage and specify which lateral section corresponds to which articular cartilage zones.Subsequently, utilization of collimated transmittance measurement and Beer-Lambert's law allowed us to estimate the μ s of articular cartilage zones.Leveraging this information alongside the optical properties of articular cartilage bulk tissue, we further calculated the g of each cartilage zone.
Furthermore, we utilized μ a and μ s of bulk articular cartilage, μ s and g of articular cartilage zones, and volume percentage of articular cartilage zones, estimated from the PLM imaging, to calculate P D avg , P D max , LS max , and P L avg of photons in articular cartilage matrix.To do so, we employed MC simulations under various scenarios of articular cartilage geometry, including single-layer and multi-layer configurations.These simulation approaches facilitated the estimation of P D avg , P D max , LS max , and P L avg of photons within the tissue, providing a comprehensive analysis of light propagation considering the optical properties of both bulk and zonal layers of articular cartilage.Collectively, this methodological strategy aimed to holistically investigate the intricacies of light propagation in articular cartilage, offering valuable insights into instrumentation requirements for developing an effective noninvasive optical device that can probe the tissue volume effectively.The technical details of each step of our methodology strategy are provided below.

A. Sample Preparation
In this study, we obtained osteochondral (cartilage on bone) samples (cylindrical plugs, n = 22) from different anatomical sites of seven bovine knee joints.These osteochondral plugs were harvested from two adjacent locations on the lateral and medial compartments of the femur, tibia, and patella of the joints.As a result, we obtained two sets of osteochondral samples, Set A (n = 11) and Set B (n = 11).The samples were extracted using a 15 mm-diameter punch and then subjected to optical spectroscopic and imaging measurements.To ensure consistency, the bone end of each plug was filed until it was parallel with the articular surface.
The joints were collected within one week of slaughter and stored intact at 4 °C until sample extraction (which was within one week of obtaining the joints).Furthermore, since the joints were acquired from a local abattoir, no ethical permission was required.During the collection process, no confounding factors, such as age and sex, were controlled.Throughout the extraction, the cartilage surface of the plugs was continuously rinsed with phosphate-buffered saline solution (PBS, pH 7.4, containing inhibitors) to prevent deterioration due to dehydration.Table I provides information about the number of samples collected in each set, their anatomical origin, the specific measurement they underwent, and the target parameters obtained from them.

B. MicroCT Imaging
To estimate the axial thickness and lateral diameter of the bulk articular cartilage samples, Set A was subjected to microCT imaging immediately after extraction.The imaging was performed using a microCT scanner (XT H 225, Nikon Metrology, Leuven, Belgium) with an initial voxel size set to 40 × 40 × 40 μm 3 .During subsequent 3D image reconstruction, the voxel size was changed to 50 × 50 × 50 μm 3 .From the reconstructed microCT images, the axial thickness and lateral diameter of the articular cartilage part of the osteochondral plugs were estimated.These physical properties are later used in the MC simulation to create a digital model of the articular cartilage tissue, a voxelized cylindrical grid with similar physical dimensions.After microCT imaging, the articular cartilage segment of the osteochondral samples was detached from the subchondral bone using a scalpel, and the samples were then preserved in PBS at −20 °C until required for optical measurements.Fig. 1(c)-(e) illustrate the physical properties from processing of the microCT images of the osteochondral samples.

C. Cryosectioning
To estimate μ t , μ s , and g of each articular cartilage zonal layers, lateral tissue sections of cartilage from depths representing Fig. 1.Panel (a) depicts an example of osteochondral sample harvesting into two sets of A and B from bovine knees, followed by a schematic of articular cartilage zonal architecture.In the superficial zone, SZ, the collagen fibers are parallel to the articular surface and the chondrocytes have an oval shape.In the middle zone, MZ, the collagen fibres are randomly aligned and the chondrocytes display a round shape.In the deep zone, DZ, the collagen fibres are perpendicularly aligned and the chondrocytes form clusters.Panels (b)-(f) show the distribution of the physical properties (thickness (mm), surface diameter (mm), and volume (mm 3 ) of the articular cartilage bulk samples, in addition to reconstructed images of cartilage from the microCT images.The voxel size of the reconstructed microCT images is 50 × 50 × 50 µm 3 .Panels (g) & (h) depict the schematics and photo of the integrating sphere setup used for the estimation of broadband µ a and µ s of articular cartilage bulk tissue.Panels (i) & (j) illustrate the schematics and photo of the collimated transmittance setup used for estimation of broadband µ t of articular cartilage zones.Panel (k) displays the execution path of the methodology employed in this study.
each zone were extracted for collimated transmittance measurement.To achieve this, samples in Set B were cryosectioned.The procedure involved mounting and fixing the bulk samples (from the bone end) on the rotary microtome with a cryostat chamber (Leica CM3050 S, Leica Biosystems, Wetzlar, Germany) using a fixative solution (Optimal Cutting Temperature (OCT) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
compound, Thermo Fisher Scientific Ltd., Runcorn, U.K.).The samples were allowed to freeze completely by maintaining a stationary state at −18 °C for 5 minutes.
The sectioning process was initiated by cutting and discarding the first 50 μm of the articular surface to provide a uniformly smooth and flat surface for subsequent sectioning.A series of cartilage layers, each with a thickness of 100 μm, was then consequently sectioned.The sections were hydrated with PBS and sandwiched between a glass slide (thickness = 1 mm, Menzel-Gläser Frosted Microscope Slides, ThermoFisher Scientific Oy, Finland) and a coverslip (thickness = 0.13 mm, Menzel Microscope Coverslips, ThermoFisher Scientific Oy, Finland).The prepared sections were then stored in a humid box at −20 °C.
For collimated transmittance measurements, specific sections from depths known to be consistent with each zone were selected to accurately represent the zonal variation of articular cartilage accurately.The selection criteria were as follows: I) the first section was always designated as SZ; II) the subsequent section was considered to represent MZ; and III) from the remaining sections, the one with the largest lateral diameter and uniform thickness was chosen to represent DZ.

D. Optical Setups for Optical Properties Estimation
The optical experiments conducted in this study comprise two stages.In the first stage, the samples in Set A were subjected to integrating sphere measurements to estimate μ a and μ s of articular cartilage bulk tissue.In the second stage, the thin cartilage sections in Set B were subjected to collimated transmittance to obtain μ t , μ s , and g of articular cartilage zones.A detailed description of the integrating sphere and collimated transmittance setups is given elsewhere [9], [17], [18].In summary, the integrating sphere apparatus is a dual-beam spectrophotometer operating in the spectral region of 400-1400 nm (Fig. 1(g)).A single integrating sphere with a reference beam is a common setup in optical measurements for estimating reflectance and transmittance of a sample.The integrating sphere is a spherical cavity coated with a highly reflective material (such as Spectralon), and it is designed to collect and diffuse light incident on its inner surface.The sphere serves to integrate the light over all angles, making it especially useful for measurements of diffuse reflection and transmission.In this setup, the sample is typically placed at the entrance port of the integrating sphere.Light from a source is directed onto the sample, and the integrating sphere collects both the reflected and transmitted light.The reference beam, on the other hand, is directed into the sphere but bypasses the sample, providing a baseline measurement without interaction with the sample.For reflectance measurements, the ratio of the reflected light collected by the integrating sphere to the reference beam gives the reflectance of the sample.Similarly, for transmittance measurements, the ratio of the transmitted light collected by the integrating sphere to the reference beam provides the transmittance of the sample.The computational algorithm for estimating μ a and μ s is an iterative optimization algorithm based on a 3D MC solver and predefined lookup tables of μ a and μ s .In the computational algorithm, for precise estimation of the light distribution, radially and axially stacked cylindrically shaped simulation volumes were considered to estimate potential side losses of the sample and to consider the impact of the cylindrical cuvette holder [19].
The computational algorithm implements the Henyey-Greenstein scattering phase function to simulate the angular distribution of scattered photons and considers the geometry of the samples to be a perfect cylinder.For this study, fixed values of 1.358 and 0.9 from the literature were used for n and g [20] as input to the computational algorithm.In addition, the estimated lateral diameter and axial thickness of articular cartilage samples from the Set A were used to build the cylindrical geometry of the samples in the algorithm.
The collimated transmittance setup consists of a broadband light source with an illumination fiber with 100-μm core diameter, a broadband spectrometer with a detection fiber with a 200-μm core diameter, two lenses with regular achromats of f-60 mm and a diameter of 25.4 mm, a diffuser disc, and an optical fiber with 1000-μm core diameter, connecting the diffuser disc to the spectrometer.The collimated transmittance measurements were carried out in the spectral band of 400-1400 nm.Fig. 1(b) shows the schematics and images of the collimated transmittance setup.
The broadband extinction coefficient, μ t (mm −1 ), of the thin cartilage sections in the Set B were estimated by applying Beer-Lambert's Law [2] on the collimated transmittance dataset as follows: Where T ( = T sample − T Dark ) is the transmittance signal of the samples subtracted by the dark current noise of the system (T Dark ) and normalized by the reference transmission signal (T 0 = T Ref erence − T Dark ).T Ref erence is the transmittance measurement when no sample is placed in the system.Lastly, d (mm) is the thickness of the samples.Given the values of μ a and μ s of bulk articular cartilage tissue, estimated with the integrating sphere setup, and the values of μ t of articular cartilage zones, the spectral values of μ s and g of articular cartilage zones can be estimated using the following equations: g is a parameter that characterizes the angular distribution of scattered light in a medium.It is a dimensionless quantity (g ∈ [−1, 1]) that describes the deviation of the scattering pattern from isotropic (uniform in all directions) scattering.A value of g = 0 indicates isotropic scattering, where photons are scattered equally in all directions, while g = 1(−1) implies perfectly forward (backward) scattering, where photons are predominantly scattered in the forward (backward) direction [3].In MC simulations, g is a key parameter that influences the angular distribution of scattered photons.It is particularly an essential factor, as it provides insights into the behavior of light within biological tissues, influencing the penetration depth Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.and spatial distribution of light.During the integrating sphere measurements, we first assumed a fixed value of g to facilitate estimation μ a and μ s of bulk articular cartilage tissue, as the computational algorithm requires this information, and no prior information was available.Since articular cartilage is a fibrous tissue, it is a common understanding that light propagation in fibrous tissues is mostly forward oriented [21].To this end we assumed g = 0.9 during the optical properties measurement which is a valid assumption [3].However, after estimating μ s of bulk articular cartilage and μ s of articular cartilage zones, we utilized the relationship between the scattering properties of the bulk tissue and its zones through Eq. ( 3).In this context, when we estimate the values of g of articular cartilage zones, we assume μ s/ μ s represents the contribution of individual cartilage zones to the scattering properties of bulk tissue.In other words, we assume μ s of articular cartilage zones is as same as μ s of bulk articular cartilage.
From the integrating sphere measurements, μ a and μ s of bulk articular cartilage samples were measured.While from the collimated transmittance measurements, μ t , and subsequently, μ s and g of thin cartilage sections were measured.

E. Polarized Light Microscopy
After integrating sphere measurement, the samples in the Set A were subjected to PLM imaging [22], [23], [24] to estimate the axial profile of collagen fibers orientation of the articular cartilage samples.The process began with placing the samples into a fixative solution containing formaldehyde (4%, Merck, Darmstadt, Germany) and ethylenediaminetetraacetic acid (EDTA, 10%, Merck, Darmstadt, Germany).Following fixation, the specimens were dehydrated and embedded in paraffin for axial sectioning.The embedded specimens were halved from the sagittal plane, and one set of three unstained deparaffinized 5-μm axial sections was obtained from each half.The axial sections were then subjected to polarized light microscopic imaging.
The PLM imaging apparatus consisted of a microscope body (Leitz Ortholux II POL, Leitz Wetzlar, Germany), a monochromatic light source (λ = 630 ± 30 nm, Edmund Optics Inc., Barrington, NJ, USA), crossed polarizers (Techspec optics XP42-200, Edmund Optics, Barrington, NJ, USA), and a monochromatic camera (pixel size 3.5 μm, BFS-U3-88S6M-C FLIR Blackfly S, FLIR Systems Inc., USA) with a 10.0 x magnification objective.The setup consists of a sample between a polarizer and an analyzer at a 90 • angle.Six PLM measurements were carried out for each sample and the results were averaged to provide a depth-wise orientation of the collagen fibers for the samples.The measurements were conducted at 21 different orientation angles in the band of [0 • , 180 • ] with a 9 • step size.A detailed description of the algorithm for estimating the collagen fiber orientation is reported elsewhere [22], [23].
Furthermore, the depth-wise angular profile of collagen fibers was subsequently processed to estimate the thickness of each articular cartilage zone according to the orientation of collagen fibers as follows: I) SZ with the orientation of collagen fibers in the range of II) MZ with the orientation of collagen fibers in the range of 30 • − 60 • .III) DZ with the orientation of collagen fibers in the range of ≥60 • .The output of PLM imaging is the angular profile of collagen fibers from the articular surface of the bulk cartilage samples to their subchondral bone interface.The PLM images were obtained from the thin sections obtained from histology procedure.The angular profile of collagen fibers allows us to understand how much of the tissue volume is occupied by different cartilage zones (SZ, MZ, and DZ) according to the classification provided above.

F. Estimation of Light Propagation Parameters
In the field of tissue optics, MC simulation tools provide a robust framework to predict how light interacts with complex biological structures, offering valuable insights into the impact of tissue optical properties, light distribution and energy deposition within the tissue, and light propagation characteristics such as photon penetration depth and path length.These simulations take into account factors such as tissue composition, scattering, and absorption coefficients.They are particularly instrumental in studying the spatial distribution of light, estimating energy deposition, and determining penetration depths and path lengths of photons.Recent advancements in MC simulations have significantly contributed to our comprehension of light-tissue interactions, fostering improvements in optical diagnostic and therapeutic applications [1].
To estimate P D avg , P D max , LS max , and P L avg of photons in articular cartilage over the 400-1400 nm spectral band, a series of MC simulations [11], [12], [13] were carried out using the estimated optical properties.These parameters were defined as follows, and their mean value for each wavelength in the 400-1400 nm range was estimated via launching and tracking the trajectories of 100000 photons in the MC simulation: where λ is the wavelength of investigation, N is the total number of launched photons, z max, i is the maximum z-position of the ith launched photon inside the tissue, x max, i is the maximum x-position of the ith launched photon inside the tissue, y max, i is the maximum y-position of the ith launched photon inside the tissue, P L j,i is the step size of the ith photon at step j before the photon gets absorbed or scattered out of the tissue, Z m,i is the z-position of the ith launched photon at the step m where Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the weight of the ith launched photon (W j,i ) drops to 63.21% of its initial weight [4], [5].
The MC simulations were carried out in the PyXopto simulation package [11], [12], [13].Articular cartilage samples were represented in the software as a 3D cylinder with axial thickness and lateral diameter equivalent to the values estimated with microCT imaging.In addition, the modified Henyey-Greenstein phase function [25] was used to estimate the angular distribution of scattered photons.The modified Henyey-Greenstein phase function is defined as follows: where pf mHG is the modified Henyey-Greenstein phase function, c is the normalized contribution of the Rayleigh scatterer to the scattering profile of the photons, θ is the scattering angle, g is the scattering anisotropy factor, and pf HG is the Henyey-Greenstein phase function.Parameter c is obtained by fitting the experimental values of μ s and μ s to the Mie-Rayleigh power law as below: where α (= μ s (λ 0 ), mm −1 ) is related to the density of the scatterers in the tissue, λ 0 (= 500 nm) is a reference wavelength used for nondimensionalization, and c is the normalized contribution of the Rayleigh scatterers to μ s .Similarly, α * ( = μ s (λ 0 ), mm −1 ), c * , and b * are the scatterer density parameter, the normalized contribution of Rayleigh scatterers to μ s , and the size parameter of Mie scatterers, respectively [17].
In addition to the choice of the scattering phase function, the MC simulation requires μ a , μ s , n, and g as input parameters.Due to the availability of the spectral values for μ s , μ s , and g and literature values of n for each zone, various scenarios representing different combinations of the optical properties were considered as input parameters for the MC simulation of light propagation in bulk articular cartilage tissue.These scenarios are outlined in Table II.In scenario 1, we simulate the optical response of bulk articular cartilage tissue.We estimate light propagation parameters for a single-layer bulk tissue with its broadband μ a and μ s and fixed values of g and n from the literature [20], [26].In scenario 2, we considered articular cartilage as a stack of three layers.Each layer represents a different zone of articular cartilage with thickness estimated from the angular profile of collagen fibers from PLM imaging.Each layer has its broadband μ s , estimated from collimated transmittance measurement and Equation (10), and μ a of articular cartilage.In addition, a fixed value of n [20] was used for all the layers, and the zonal average of g -estimated by averaging the broadband g of articular cartilage zones were used for all three layers.In scenario 3, similar to scenario 2, we considered articular cartilage as a stack of three layers representing its zonal layers.For each zone, we assumed its μ a to be equivalent to bulk tissue μ a , and its μ s and g corresponding to the μ s and g of articular cartilage zonal layers.Lastly, a fixed value of n for each articular cartilage zone was obtained from the literature and used [20].
In the MC simulation, the measurement geometry, including the sample holder, was not implemented.This is because the light propagation inside of the tissue is of importance and mimicking the measurement geometry does not change this phenomenon.Moreover, reflectance and transmittance were simulated over the surface of the samples with their surface diameter and thickness measured by microCT imaging to match the samples used in the integrating sphere measurement as closely as possible.
Furthermore, the volume of the samples was considered as a 3D voxel mesh, with each voxel having a dimension of 50 × 50 × 50 μm 3 .In the scenarios with articular cartilage as a 3D layer model, the voxelized mesh was divided into three horizontal layers with a thickness equivalent to the thickness of articular cartilage zones estimated via PLM imaging.Fig. 2 depicts the optical properties used in each simulation scenario.

G. Statistical Analysis
To investigate whether different scenarios of the Monte Carlo simulation yield statistically different results in terms of the parameters of interest, P D avg , P D max , LS max , and P L avg , as well as simulated reflectance and transmittance, a series of statistical tests were conducted.The statistical test comprises of Kolmogorov-Smirnov normality test, Leven's test for equality of variances, One-way ANOVA and Tukey's tests to examine the difference in normally distributed observations, and Kruskal-Wallis and Dunn's tests for nonparametric observations.All computational and statistical analyses performed in the present study were carried out in MATLAB (R2019b and R2020b) and Python v3.7.

III. RESULTS
The P D avg , P D max , LS max , and P L avg over all simulation scenarios are illustrated in Fig. 3.Both P D avg and P D max exhibit an increasing trend, whereas P L avg shows a decreasing trend.P D avg was generally lower than P D max across the spectral band of 400-1400 nm.In all simulation scenarios, P D avg and P D max increased at a more significant rate in the visible range (400-700 nm) than in the NIR range (700-1400 nm).In the NIR range, the increase in the values of P D avg and P D max declined with apparent plateauing of P D avg over the NIR band.Spectral trend of LS max showed an increasing behavior over the spectral band of 400-900 nm and decreasing features over the spectral band of 900-1400 nm.On the other hand, the broadband values of P L avg decrease at a slower rate in the visible than in the NIR range.
In scenario 1, we compared the simulated measured reflectance and transmittance values and observed a relative difference of 37.79% and 14.5% for reflectance and transmittance, respectively (Fig. 4(a) and (b)).These results serve as a baseline to assess how accounting for the zonal architecture of articular cartilage can change the measured optical response of the tissue.
In scenario 2, a fixed value of n [20] was used for all the layers, and the same broadband values of g -estimated by averaging the broadband g of articular cartilage zones were used for all three layers.The simulated reflectance and transmittance exhibited a relative difference of 40.53% and 17.04% from their experimental counterparts, respectively (Fig. 4(c) and (d).When the simulated optical response of articular cartilage and P D avg , P D max , LS max , and P L avg of photons were compared to their counterparts from scenario 1, and the statistical tests yielded no significant difference between these values.
In scenario 3, the simulated reflectance and transmittance exhibited a relative difference of 38.45% and 15.48% from the experimental reflectance and transmittance, respectively (Fig. 4(e) and (f)).However, the statistical tests did not show any meaningful difference between the spectral values of the P D avg , P D max , LS max , and P L avg of photons estimated in this scenario and those estimated in the other simulation scenarios (Table III).
To test the validity and fidelity of the simulation results, the simulated reflectance and transmittance values were compared against their experimental counterparts (Fig. 4).In addition, a series of statistical tests were performed to investigate whether the simulation scenarios lead to different values of reflectance, transmittance, P D avg , P D max , LS max , and P L avg for each Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.II.Values are presented as mean ± standard deviation.Fig. 5. Visualization of photon propagation and absorption in articular cartilage for an arbitrary sample at 800 nm.The image depicts the photon weight deposition on the x-z plane from the surface of articular cartilage to its bottom end.The energy of photons in MC simulation is their initial weight deposited in each voxel of the simulation grid according to the µ a of tissue at the wavelength of investigation.The MC simulation was done according to the specifications of scenario 1, outlined in Table II.
wavelength of 400-1400 nm, when different simulation scenarios were compared with each other.Consequently, no meaningful statistical difference was found, and a summary of the statistical test is presented in Table III.Additionally, a visualization of photon propagation and their absorbed energy in articular cartilage bulk tissue (scenario 1) are provided in Fig. 5 to enhance the understanding of how photons are scattered and absorbed in the matrix of articular cartilage.When thickness of bulk tissue and zones of articular cartilage samples used in this study is considered, on average, bulk tissue has a thickness of 1.4857 mm, SZ has a thickness of 0.1128 mm, MZ has a thickness of 0.2683 mm, and DZ has a thickness of 1.1046 mm.When these values are compared with the P D avg and P D max of photons in the visible and NIR bands, the results indicate that visible and NIR photons travel beyond the SZ and MZ layers of articular cartilage and propagate into DZ.However, the photons in the visible band see less of articular cartilage DZ than the NIR photons (Fig. 3(a)-(f)).
Furthermore, when LS max and P L avg of photons are considered, photons in the visible band show a high amount of distance travelled, yet small lateral spread and penetration depth.This result indicates that photons in the visible range (400-900 nm) probe a localized superficial volume of the articular cartilage tissue that are mostly limited to the SZ and MZ layers of the tissue.The peak values of LS max are observed in the spectral region of 700-1100 nm, combined with a penetration depth of 1 mm, it seems that this spectral region can probe the full matrix of articular cartilage.On the other hand, LS max values show a decline over the spectral region of 1100-1400 nm.However, the penetration depth in this region plateau to a value around 1.0-1.2mm, suggesting the decrease in P L avg is due to the decrease in LS max of photons in this band (Fig. 3(g)-(l)).
Around 1400 nm, a sudden drop (knee point) in the values of P D avg , P D max , LS max , and P L avg is observed.This is mainly due to the strong absorption of light by articular cartilage tissue in this spectral band.Fig. 2(a), (c), and (f) suggest a 10-fold increase in μ a of articular cartilage, compared to the spectral band of 1200-1300 nm.Water is the main chromophore in articular cartilage in the NIR spectral region.Therefore, it is expected that the sudden decrease in the values of P D avg , P D max , LS max , and P L avg is because of strong light absorption by the interstitial fluid of the tissue.

IV. DISCUSSION
The knowledge of P D avg , P D max , LS max , and P L avg of photons in tissues, in the context of optical-based medical devices for orthopedic diagnostics, especially in articular cartilage, is critical for optimizing device design and enhancing diagnostic capabilities.Articular cartilage, comprising distinct zones with varying composition and structure, presents a unique challenge for instrumentation due to its layered architecture.Understanding penetration depth, lateral spread, and pathlength of photons is crucial for the following reasons: I) Zone-specific information retrieval: articular cartilage consists of different zones, each with specific structural and compositional features.The ability of optical devices to penetrate into these zones depends on the penetration depth.This information allows for the design of instruments that can effectively retrieve information from different layers, providing a more comprehensive understanding of the cartilage structure and composition.II) Optimization for diagnostic accuracy: the optimization of device design based on penetration depth, lateral spread, and pathlength information enables improved diagnostic accuracy.Devices can be tuned to specific wavelengths to capture relevant data from specific depths within the cartilage, ensuring that the diagnostic information obtained is representative of the targeted region.This is crucial for detecting and monitoring degenerative diseases like osteoarthritis, where disruptions in tissue homeostasis can vary across different zones.III) Tracking disease onset and progression: Osteoarthritis often manifests with changes in specific regions of articular cartilage.Devices optimized for varying penetration depths can track the onset and progression of degenerative diseases by providing detailed information about changes in tissue composition, structure, and optical properties at different depths.This is essential for early diagnosis, personalized treatment planning, and monitoring disease progression over time.
In this study, we reported, for the first time, the P D avg , P D max , LS max , and P L avg of photons in articular cartilage in the visible and NIR spectral regions (400-1400 nm).We achieved this by utilizing MC simulation [11], [12], [13] and optical properties of articular cartilage zonal layers and bulk tissue [9], [10], [20].Our results indicate that photons in different regions of the visible-NIR light behave differently.We demonstrated that the photons in the 400-700 nm penetrate only to the SZ and MZ layers of the tissue with a low lateral spread, suggesting localized probation of tissue volume.However, in the 700-1100 nm we observed a maximal lateral spread of photons and approximate penetration depth of 1 mm.This means that photons in this band, can reach the DZ layer of articular cartilage and spread as wide as 6 mm in diameter which is optimal for wide and deep probation of tissue volume.Lastly, the photons in the 1100-1400 nm ban, showed a maximal penetration depth while their lateral spread declined compared to the 700-1100 nm, suggesting a narrow yet deep probation of tissue volume that can reach beyond the SZ and MZ layers of articular cartilage and probes the DZ layer of the tissue.These findings support the applicability of diffuse optical modalities for probing the volume of articular cartilage and providing indicators of its health.
In our analysis, we investigated how the optical response and characteristics of light propagation in articular cartilage vary according to three scenarios (Table II).Our goal was to explore the effect of optical properties, especially μ s and g, of articular cartilage zones, on the reflectance and transmittance of its bulk tissue, as well as the P D avg , P D max , LS max , and P L avg of photons.Our results indicate that a single-layer representation of articular cartilage tissue can potentially describe light propagation in its matrix in the visible range (Fig. 4(a) and (b)).A three-layer model of articular cartilage, based on the scattering properties of its zonal layers, provided similar accuracy in describing the light propagation in the layered matrix of articular cartilage (Fig. 4(c)-(f)).
When the P D avg , P D max , LS max , and P L avg of photons in articular cartilage in the visible are considered (Fig. 3), our results suggest that the visible photons can reach only the SZ and MZ layers of articular cartilage (Fig. 3(a), (d), and (g)).In contrast, the photons in the 700-1100 nm band can reach deep into the tissue and probe its full depth, while simultaneously exhibit a relatively large lateral spread.As a result, in diagnostic applications, when measuring the full depth of articular cartilage is required, the photons in this band can provide a better assessment than visible photons (Fig. 3(b), (e), and (h)).Furthermore, when assessment of a localized and specified volume of a superficial region of articular cartilage is required (for example around a surface defect), visible light can potentially provide a more comprehensive assessment due to their localized trajectory in the superficial part of the tissue (Fig. 3(c), (f), and (i)).
In this study, certain reasonable assumptions were made.We assumed that light propagation in articular cartilage follows an isotropic light propagation model and that the optical properties and response of the tissue are independent of the angle of incident light.The underlying reason is due to numerous unknowns associated with the anisotropic model of light propagation in articular cartilage such as the microscale orientation of fibers Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.and change in their diameter from the surface to the bottom of the tissue, that makes it non-feasible [1].The second assumption was that the μ s and g of articular cartilage zones are applicable for estimating the optical response of bulk articular cartilage tissue when a layered model of the tissue is studied.In theory, the three-layer model of articular cartilage should have provided a better accuracy in describing the light propagation in the tissue, compared to the one-layer model.However, due to intricacies related to the orientation of the collagen fibers, similar μ a for all the layers, and negligence of the impact of tissue cell population, the assumption of tissue being stack of three slabs did not improve the results obtained when tissue was considered a single layer.However, the three-layered model provided same level of accuracy as the single-layer model.
Furthermore, when possible sources of uncertainty and inaccuracy in this study are considered, they are accumulation of inaccuracy associated with each step of the employed methodology, including sample harvesting, microCT imaging, integrating sphere measurement, collimated transmittance measurement, and polarized light microscopy imaging, and MC simulation.Kafian-Attari et al. [10], [26] have already addressed the degrees of uncertainty associated with microCT imaging technique, integrating sphere measurement, and collimated transmittance measurement utilized in this study.Additionally, the polarized light microscopy imaging and the protocol we utilized for sample preparation and imaging are based on the standard protocols disclosed and outlined in the literature [23].To estimate the degree of uncertainty associated with MC simulation, we repeatedly performed the MC simulation for one arbitrary sample at 410 nm five times and recorded the simulated values of reflectance and transmittance.The mean ± standard deviations of the simulated values of reflectance and transmittance were 78.66% ± 0.1436% and 6.483 ± 0.3082%, respectively.Thus, it seems that for large values of reflectance and transmittance, the deviation observed in the results is small.However, as the reflectance and transmittance tend to become smaller, the impact of deviation due to the MC simulation becomes apparent and pronounced.Given that this is the first time such a study has been conducted in this domain, there might be some unidentified biases and sources of uncertainty that affect the results presented in this study.Yet we firmly believe that this study is a forward step in the right direction of understanding the importance of light propagation parameters for instrumenting optimized optical device that can positively impact the diagnosis-treatment pathways in orthopedic surgeries.

V. CONCLUSION
This study reports the broadband light propagation characteristics, P D avg , P D max , LS max , and P L avg , in articular cartilage in the visible-NIR region (400-1400 nm).Our findings suggest that the photons in the visible band probe the lateral volume of articular cartilage SZ and MZ while the photons in the NIR band probe the deeper tissue.Moreover, we showed that a single-layer model of articular cartilage can reliably describe the light propagation in articular cartilage.This indicates that diffuse optical spectroscopy can potentially be adapted for objective assessment of articular cartilage tissue.

DISCLOSURE
The Authors declare no conflict of interest.

Fig. 2 .
Fig. 2. Input optical parameters for simulation scenarios: the broadband µ a and µ s of bulk articular cartilage used in scenario 1 (a & b); the broadband µ a , zonal µ s , and the broadband values of g, averaged over the zones, used in scenario 2 (c-e); and the broadband µ a , zonal µ s and g used in scenario 3 (f-h).SZ, MZ, and DZ are superficial, middle, and deep zones of articular cartilage, respectively.µ a (mm −1 ), µ s (mm −1 ), and g are the absorption coefficient, scattering coefficient, and scattering anisotropy factor, respectively.Values are presented as mean ± standard deviation.

Fig. 3 .
Fig. 3. Panels (a)-(c) depict the average penetration depth (P D avg , mm) of photons over simulation scenarios 1-3.Panels (d)-(f) show the averaged maximum penetration depth (P D max , mm) of photons over simulation scenarios 1-3.Panels (g)-(i) illustrate the average path length (P L avg , mm) of photons over the spectral band of 400-1400 nm over simulation scenarios 1-3.Panels (j)-(l) show the averaged maximum lateral spread of photons (LS max , mm) three simulation scenarios outlined in Table II.Values are presented as mean ± standard deviation.

Fig. 4 .
Fig. 4. Comparison of simulated and experimental reflectance (R, %) and transmittance (T , %): according to the simulation scenario 1 (a) & (b), scenario 2 (c) & (d), and scenario 3 (e) & (f).The mean relative difference between the simulated and measured values of reflectance and transmittance is depicted with a red curve.The details of the simulation scenarios are outlined in Table II.Values are presented as mean ± standard deviation.

TABLE I DETAILS
OF THE SAMPLES COLLECTED IN SETS A AND B

TABLE II DETAILS
OF THE INPUT PARAMETERS USED IN THE SIMULATIONS

TABLE III THE
P-VALUE OF POSTHOC METHOD OF GROUP STATISTICAL TEST FOR ASSESSING THE STATISTICALLY SIGNIFIANT DIFFERENCE BETWEEN THE SIMULTION RESULTS