1. Introduction

Traditional techniques like arthroscopy, clinical magnetic resonance imaging (MRI), radiography and polarized light microscopy as diagnostic tools to study the onset and development of osteoarthritis have not been completely successful and each has its own drawbacks [1–3]. A non-invasive imaging tool to obtain information about the microstructure of the collagen fiber organization in articular cartilage has potential application in understanding and diagnosing different phases of osteoarthritis. Optical coherence tomography (OCT) is a promising candidate in this regard with its non-invasive ability and good resolution over penetration depths of 1-2mm [4]. Specialized versions of OCT which make use of polarized light: polarization–sensitive optical coherence tomography (PS-OCT) offer more information on the structural orientation of the collagen fibers in articular cartilage [5,6]. Osteoarthritis is associated with the disruption of the delicate 3D structural organization of the Type II collagen fibers in the hyaline cartilage. Collagen fibers in articular cartilage tissues are organized in a complex 3D structural organization. A detailed study of the birefringence information obtained from light beam incident on the sample at multiple angles of incidence provides a more complete picture of the 3D orientation of the collagen fiber fast axis. This has been first shown by Ugryumova et al., for a single layered model of collagen fibers in tendon tissues [7]. A more rigorous theoretical basis to study the propagation of the polarized components of light in layered anisotropic media is based on the extended Jones matrix calculus (EJMC), developed by Yeh et al., for describing light propagation in liquid crystals. The EJMC describes light propagation through layered birefringent media with an arbitrary 3-D orientation of the optic axis in each layer by extending the conventional 2 × 2 Jones calculus to the case of off axis light transmission, accounting also for the Fresnel transmission coefficients at the initial interface [8,9]. This provides a framework for studying layered polarized light propagation in articular cartilage, which possesses a complex 3D orientation of collagen fibers. The EJMC has been used by our group previously to theoretically model light transport in a scattering medium using a Monte Carlo method [10]. Also, a purely theoretical description pertaining to cartilage using the EJMC has been reported recently by Fanjul-Velez et al. [11]. They have also applied this approach to carry out a polarimetric study of the human cornea over its layered structure.

In this paper, angle- and depth-resolved retardance profiles obtained from bovine articular cartilage using a time-domain PS-OCT system are interpreted using a layered model, which is based on scanning electron microscopy (SEM) data [12] and then propagated through the EJMC. The 3-D geometry of the collagen fibers is then crudely estimated by parameterizing it and then manually adjusting these parameters to obtain the least-squares best fit of the data. To the best of our knowledge, this is for the first time a comparative study based on experimental and theoretical data of microstructure of articular cartilage has been reported in order to explain the PS-OCT data obtained.

2. Materials and Methods

The layered morphology of the collagen fiber network in articular cartilage is as shown in Fig. 1(a) and also shown in more detail in Fig. 2(b) . The surface ‘tangential’ layer comprises typically 10% of the total thickness of the articular cartilage and consists of collagen fibers oriented parallel to the surface. This layer is followed by the ‘transitional’ zone, which as per the lamellar model of Jeffery et al. has collagen fibers organized into parallel sheets which arch downwards leading to the radial zone, in which the collagen fibers are oriented perpendicular to the surface [12]. The transitional zone comprises 40-60% of the total thickness of the cartilage with the rest of the thickness comprising of the radial zone leading to the subchondral bone [13]. We base our theoretical model of the articular cartilage on this leaf-like lamellar model of Jeffery et al., starting with a simple model with the assumption of constant azimuth angle orientation of the collagen fibers throughout the depth of the tissue [12]. Fresh tissue samples of bovine articular cartilage were extracted from the fetlock joint of the hindlimb of the animal obtained from the local abattoir. The samples were stored frozen at −20°C prior to imaging and then sectioned along the anterior side of the apex for imaging. Multi-angle (0°, +60° and −60° relative to the surface normal) PS-OCT measurements were carried out for two different orthogonal planes, $x-z$and $y-z$ (Fig. 1(b)). The $x,y,z$Cartesian coordinate system is defined such that the z-axis represents the axial depth direction. The $x-z$and$y-z$planes corresponds to the coronal and sagittal planes, defined as the plane in which the two connected bones are constrained to move and the orthogonal plane, respectively. Hence in general, these planes are oriented at an unknown angle with respect to the azimuthal orientation of the superficial collagen fibers.

 

Fig. 1 (a) A schematic showing the EJMC approach to modeling the zonal layered structure of articular cartilage. Here ϕ_c and θ_c give the azimuthal and the polar angle of the collagen fiber direction, i.e., the optic axis. The incident angle θ_inc is varied with respect to the surface normal of the tissue over two orthogonal planes to obtain the multi-angle PS-OCT data. 1(b) A schematic of the imaging procedure of multi-angle PS-OCT study of articular cartilage over two orthogonal planes.

Download Full Size | PDF

 

Fig. 2 (a) A schematic of the bulk-optics TD-PS-OCT system used to carry out the study; QWP, quarter wave plate, BS, beam splitter, PBS, polarizing beam splitter, P, polarizer. (b) A schematic of the cartilage zonal layered structure and the layer thickness used for the EJMC study. Also shown are the orientations of the polar angle of the collagen fast axis varying from 90° in the superficial zone to gradually becoming 0° in the radial zone.

Download Full Size | PDF

The theory of the time-domain PS-OCT (TD-PS-OCT) system used for our study is based, is similar to the first PS-OCT system demonstrated by Hee et al. [5], and other first generation TD-PS-OCT systems reported subsequently [14]. Our system is described in details elsewhere [15]. A schematic of the experimental set up is as shown in Fig. 2(a). The bulk optics TD-PS-OCT system uses a light source based on quantum dot superluminescent diode centered at 1300nm with a bandwidth of 85nm [16]. The light beam is linearly polarized and then split into the reference arm and the sample arm. In the reference arm, a quarter wave plate is placed with fast axis oriented at 22.5°, which on forward and backward propagation along the reference path yields orthogonal components of equal amplitude. In the sample arm, the linearly polarized light is incident on a quarter wave plate oriented at 45° yielding light incident on the sample that is circularly polarized. The Fourier delay line scanning speed is 100Hz which means it takes 4s to acquire an OCT image of 400 lateral scans (A-scans) over a width of 4mm. From the intensities detected at the two detectors placed after a polarizing beam splitter, $I_V(z)$and$I_H(z)$, the depth dependent retardance values are obtained as

The nondepolarizing nature of the backscattering measurement involved in OCT allows Jones matrix based analysis of the measured signals [17]. Henceforth, we base our simulation on Jones matrix analysis, incorporating the extended Jones matrix calculus which takes into account off axis light propagation in birefringent multi-layered structures. We can model the light beam path of the circularly polarized light incident on the sample and then back onto the detectors through Jones calculus as

where $A_H(z)$and $A_V(z)$are the electric field horizontal and vertical amplitudes of the light beam on the detectors. $QWP$represents the quarter wave plate inserted in the sample. $J_{sample}$is the round-trip Jones matrix for the layered cartilage model, which we calculate using the EJMC. EJMC can be applied to the birefringent-layered cartilage model on the assumption of weak birefringence ($\Delta n=\left|n_e-n_o\right|<<n_o,n_e$), which is satisfied for biological tissues such as cartilage. $n_e$and $n_o$are the refractive indices of extraordinary and ordinary modes of light propagation, respectively. Here, we present the analysis based on the electric field components of the two orthogonal components of the light beam detected. Although in principle, phase-resolved data analysis can yield complex depth-resolved amplitudes, in practice, the depth dependent retardance profiles are calculated using only the modulus value of the amplitudes, which then restricts the values of the retardances calculated to the range $\left[\begin{array}{cc}0& \pi /2\end{array}\right]$for the experimental data.

With the assumption of negligible multiple reflection occurring at the interface of the different layers in the cartilage model, the overall layered cartilage model is given by

where $T_i=\left(\begin{array}{cc}{t}_s& 0\\ 0& t_p\end{array}\right)$and $T_o=\left(\begin{array}{cc}{t}^{\text{'}}{}_s& 0\\ 0& t^{\text{'}}{}_p\end{array}\right)$are the Fresnel reflection coefficients at the air tissue interface and P is the single-pass cumulative Jones matrix of the cartilage from the surface to the desired depth. For a structure comprising of m parallel layers, P can be represented as a product of Jones matrices of individual layers with each layer treated as a linear retarder with a fixed fast axis orientation:

where $R({\psi }_i)=\left(\begin{array}{cc}\cos {\psi }_i& -\sin {\psi }_i\\ \sin {\psi }_i& \cos {\psi }_i\end{array}\right)$gives the transformation matrix for the rotation of the coordinate system from that describing the normal modes of propagation (i.e., the o-wave and e-wave) to that describing the s-wave and p-wave. The full equations for $T_i,T_o,R({\psi }_i),k_{oz}$ and $k_{ez}$ can be found in [8] and [9]. In our model, the layered structure of cartilage is described by dividing the full cartilage thickness into$m=22$ layers each of thickness $d_i$(Fig. 2(b)). The superficial layer is assumed to be 100μm thick divided into 4 layers of uniform thickness. The transitional and radial layers are 500μm and 400μm respectively with each further divided into layers of thickness 50μm. The polar angle of the superficial zone is 90° which gradually changes in the transitional zonal layers becoming 0° in the radial layer. The polar angle change across the depth of the tissue is based on the SEM studies [12]. In this study, the values of the polar angle input for the transitional and radial zone (for the odd numbered layers) are [85, 60, 45, 35, 28, 15, 0, 0, 0] degrees. All layers are assigned the same azimuthal fiber orientation along the depth of the tissue, an arrangement suggested by the lamellar arrangement of the fibers shown in the SEM images. A small azimuthal angle in the range of 0-10° seems to provide a closer fit to the PS-OCT data than zero (azimuthal angle of 0° was used for this simulation); recall that the azimuth is unknown a priori as it is well known that the orientation of split-lines in articular cartilage is not necessarily parallel to the sagittal ridge. The refractive indices of the ordinary and extraordinary modes of light beam propagation are assumed to be constant over the depth of the cartilage. Whilst this assumption could be questioned, our aim in this study is a preliminary investigation to establish the applicability of this type of modeling to articular cartilage and many refinements of the model are possible due to its inherent flexibility. The values of birefringence reported in the literature varies from $(0.45\pm 0.12)\times {10}^{-3}$for 4 month old porcine sample to $(1.48\pm 0.55)\times {10}^{-3}$for a 21 month old porcine sample [18] as well as a value of $3\times {10}^{-3}$reported by Jiao et al., for bovine sample [19]. The theoretical model previously applied by Fanjul-Velez et al., assumes a high value of birefringence ($~5\times {10}^{-3}$) for articular cartilage, which is much higher than the value we find for bovine cartilage experimentally [11]. A value of $\Delta n=2.2\times {10}^{-3}$gave a good manual fit to the experimental data for the other parameters as stated. The polar angle of the fibers varies in the manner described above. All these parameters are propagated through Eqs. (2)–(4) to derive $A_H(z)$ and $A_V(z)$. From these quantities the single-pass phase retardance in the noise-free case could be calculated using Eq. (5) [5]:

In practice, however, Eq. (5) cannot fully describe the experimental data because of the effects of measurement noise. A noise model based on the OCT signal attenuation in the single scattering approximation $I(z)\propto \exp (-2{\mu }_{t}z)$is introduced. Here, ${\mu }_t$is the attenuation coefficient. A noise bias term σ is also introduced into the noise model similar to the background noise reported by Schoenenberger et al. [20]. The modified equation for the depth-dependent retardance is then given as

3. Results and Discussions

For the above mentioned choice of parameters, angle-resolved depth-dependent retardance profiles for oblique incidences of ±60° as well as for normal incidence in two orthogonal planes in the chosen coordinate system are simulated using the EJMC and closeness to the experimental data is visually assessed to obtain a good manual fit. This procedure of manual fitting was chosen initially to assess the likely success of an automated fitting algorithm based on nonlinear optimization.

The simulated profiles are compared with the experimental data in Figs. 3(a)-(e) as obtained by fitting the parameters described and visually assessed for 5measurements of angle-resolved PS-OCT measurements together. Note that the measured profiles are obtained by averaging 50 neighboring A-scans to reduce the speckle noise. The range of the retardance value is$\left[\begin{array}{cc}0& \frac{\pi }{2}\end{array}\right]$, as required by the inverse tangent function in Eq. (5) and the positivity of$A_V(z)$ and$A_H(z)$. Wrapping of values that lie outside this range creates the obtained retardance profiles. Inclusion of the noise model as given in Eq. (6) makes the retardance profiles wrap before reaching the extremes of this range since neither the numerator nor the denominator in Eq. (6) can reach zero in the presence of noise bias. As plotted in Fig. 3(c), for normal incidence as the incident beam direction is perpendicular to the collagen fiber optic axis direction in the superficial zone but parallel in the deeper radial zone, the cumulative retardance tends to display a strong band near the surface and remain constant in the deeper region, an observation which is supported by EJMC simulation. Figures 3(a)-(e) reveals how the banding pattern for retardance profiles over oblique angle of incidence along the two orthogonal planes of imaging obtained from PS-OCT experimental data tend to agree with the results obtained from EJMC simulation for our assumed cartilage model. A reasonable manual fit is obtained between the predicted retardance profiles and the experimental data. Note how the plane in which sheets are folded (i.e. the plane containing the surface normal of the sheets) can be located since incident illumination directions lying in this plane yield the maximum difference in banding periodicity (Figs. 3(a) and 3(b)).The result presented here has to be further supported based on histology studies to fully validate this technique as a tool to study the microstructure of articular cartilage. This work is ongoing in our lab.

 

Fig. 3 (a)-(e) Depth-dependent retardance profiles obtained from PS-OCT measurements (dots) compared with those obtained from EJMC simulation results (line + star) over multiple angles of illumination in sagittal and coronal planes of the bovine cartilage sample using manual fit.

Download Full Size | PDF

To extend the technique towards potential clinical applications an automated fitting procedure has also been developed using the simplex optimization algorithm ‘fminsearch’ from optimtool toolbox of Matlab^TM. It is based on Nelder-Mead direct search optimization method which works by minimizing an objective function [21]. The routine optimizes 8 free parameters: the thickness of superficial zone, $t_{\sup }$, the thickness of transitional zone, $t_{trans}$, the quadratic and linear coefficients ‘a’, ‘b’ of a quadratic function to describe the arching of the collagen fibers in the transitional zone, the true birefringence, $n_e-n_o$, the azimuthal orientation of the optic axis, ${\varphi }_c$, the tissue attenuation coefficient, ${\mu }_t$, and the noise bias term, σ. The error metric (Eq. (7)) that is minimized is

where n denotes the number of axial measurement points, and j labels the incident k-vectors. The depth dependent change of ${\theta }_c$is modeled as a quadratic profile (Eq. (8)) in the transitional zone as

where z represents the depth into the transitional layer. The initial values of a, b, and c are obtained from the tabulated ${\theta }_c$values used in the manual fit described earlier. Also, the value of ${\theta }_c$ gradually varies from 90° to radial along the transitional layer which restricts the coefficient ‘c’ to small bound [85 90] which typically gives ${\theta }_c$value at the start of the transitional layer. However, in this study, coefficient ‘c’ in the quadratic profile of ${\theta }_c$ is not a variable input into the optimizer, but is input as a constant of value 90 which describes ${\theta }_c$at the onset of the transitional zone. The list of the parameters used in the optimizer and their initial and final values are tabulated in Table 1and shown in Figs. 4(a)-(e) . A noticeable improvement in the fit is produced by the optimizer although significant differences still persist and indicate that the model is still in need of further refinement.

 

Fig. 4 (a)-(e) Depth-dependent retardance profiles obtained from PS-OCT measurements (dots) compared with those obtained from EJMC simulation results (line + star) over multiple angles of illumination in sagittal and coronal planes of the bovine cartilage sample based on parameters obtained from ‘fminsearch’ optimizer and bootstrap Monte Carlo method.

Download Full Size | PDF

The optimizer can also be used to estimate the precision on the fitted parameters, using the “Bootstrap Monte Carlo” procedure. We apply the optimizer to an ensemble of 100 random data sets, generated using the mean experimental data set and the estimated standard deviations at each axial point [22]. The statistics of the fitted parameters values are then measured directly. The mean profiles are averaged over 50 neighboring A-scans. To estimate the standard deviations we extract 10 A-scans randomly from the original set of 50 and generate one realization of the average. This is then repeated 10 times and the standard deviation at each axial pixel measured. The standard deviation for the 50 A-scan averages is then taken to be 1/√5 of this value. The resulting mean values along with the standard deviations of the 9 fitted parameters are tabulated in Table 1 . It is also noted that, within the set of constraints to the parameters given by physically possible values, a reasonable choice of any input conditions yield the output within the error range presented.

Table 1. List of all the parameters used in the ‘fminsearch’ optimizer for obtaining a good fit of the depth-dependent retardance profiles obtained from angle-resolved PS-OCT imaging and those obtained from EJMC simulation^a

View Table

In summary, we have shown the applicability of the EJMC approach to understanding the complex 3D structure possessed by articular cartilage by analyzing the PS-OCT data obtained for multiple angles of illumination. Although the PS-OCT system that we have used is not state of the art, the analysis shown in this paper is quite general. Further studies are underway to use this EJMC approach to quantitatively map the topographical variations in the 3D structural over the entire articulating surface, which ultimately could help in better understanding the biomechanics of the tissue and also detecting pathological changes. We have also applied this EJMC simulation to PS-OCT data obtained from a recently implemented continuous polarization modulation swept source PS-OCT system, the speed of which could potentially allow these measurements to be made in vivo. Results will be reported in due course.