1. Introduction

Light propagation in bulk tissue, i.e. the “photon migration” process, is usually described by the radiative transport theory, or its diffuse approximation [1] in many cases. For tissues with isotropic structures, photon scattering is independent to the incident direction within the medium. Therefore, the resulting diffuse reflectance at sample surface is rotational symmetric when point-incident light is used. However, in some tissues, such as dentin and skin, the highly organized structures produce anisotropic scattering patterns [3–8]. Several methods have been developed to study light propagation in such anisotropic media. Kienle et al applied the Mie scattering of infinite cylinder in a Monte Carlo model to study photon migration in dentin and artery tissues [3, 4]. Heino et al [5] introduced a diffuse tensor to replace the usual scalar coefficient in diffuse equation to describe the anisotropic diffusion process. Based on continuous time random walk (CTRW) theory, Gandjbakhche et al [6, 7] derived analytical solutions of light reflectance and transmittance in anisotropic media. For turbid medium with a preferred light scattering orientation, both Heino’s and Gandjbakhche’s models predicted an elliptical reflectance profile, which was observed experimentally in several anisotropic tissue models [8–10].

One extraordinary anisotropy tissue that has not been thoroughly studied in the field of tissue optics is muscle. Striated or skeletal muscle is the single most abundant tissue in most vertebrates and is essential for locomotion and other important physiological functions. The normal mechanical and physiological functions of all striated muscles are realized by highly organized structures termed sarcomeres [11]. Sarcomeres are aligned precisely in muscle fibers and are readily observed using light microscopy as alternating light and dark bands called I-band and A-band. The periodic appearance of I- and A-band is due to the difference in their optical refractive indices. Such unique structure makes laser diffraction a useful tool in studying sarcomere organizations in single fiber [12, 13]. However, it is not clear whether the effect of sarcomere structure is still significant for multiply scattered light in whole muscle.

Because of the inherent fibrous structure, light propagation in muscle is anisotropic [14], and optical properties of muscle depend on measurement orientation [15, 16]. More importantly, our recent studies [17] have found that light scattering in bulk muscle is strongly affected by sarcomere structure. In this paper, we investigated the effect of sarcomere structure on 2D optical reflectance distribution in whole pre-rigor muscles. We discovered a unique reflectance image profile that has not been reported before. In order to quantify the reflectance images, we developed a mathematical model to fit the equi-intensity profiles. The quantitative fitting parameters were used to characterize the changes in reflectance images during stretching and rigor process. None of existing models can explain our observations. However, we demonstrated that such unique image patterns can be produced when incorporating the sarcomere diffraction in a Monte Carlo model.

2. Materials and methods

2.1 Sample preparation

All muscle samples were obtained from the Meat Science Laboratory at the University of Missouri-Columbia. Sternomandibularis muscle (~5×20×5cm) was chosen in this study because it can be excised from beef cattle within 5–10 min after the slaughter. It has a uniform muscle fiber orientation, which is a good choice for our studies. All external fat layers around the muscle were carefully removed after sample acquisition.

In addition to muscle samples, isotropic scattering phantoms were made from 20% Intralipid^® (Sigma Aldrich, Inc., St Louis, MO) solution to verify the linearity of CCD and fitting algorithms. Protein extrusion products obtained from high moisture extrusion [10, 18] of soy proteins were used to represent fibrous tissues with strong anisotropic structures.

2.2 Experimental setup

Muscle samples were placed on the sample holder so that the muscle fiber direction was approximately aligned with the x-axis (Fig. 1). Each end of the muscle sample was clamped firmly during the experiment. A black plate with a glass (thickness = 0.13~0.16mm, n=1.52) window (~3.5×4.5 cm²) was placed above the sample to create a flat surface. An optic fiber (400 μm) was held slightly above the window at an oblique incident angle of ~43–47° to the sample surface. It was located within the plane formed by the y-axis and the normal vector to the sample surface. The other end of the optic fiber was coupled to a red LED (λ=680nm, 2.2V, 350mA). The incident light intensity on the sample surface was ~120 μW.

 

Fig. 1. A schematic illustration of the experimental setup. A laboratory coordinate system was set so that the x-axis is along the muscle fiber orientation. Light was incident obliquely via a 400 μm fiber.

Download Full Size | PDF

A digital camera (Model 1010, JAI-Pulnix, San Jose, CA) was mounted above the sample. The camera was aligned so that the light incident point was approximately at the center of the acquired image. Images were transferred to a personal computer via a frame grabber (Meteor-II, Matrox Inc., Quabec, Canada). In order to improve the signal to noise ratio, 10 consecutive images were recorded and averaged after subtracting the dark current. The reflectance image was stored in a 1001 × 1016 array for further processing.

When studying the effect of rigor process, muscle samples were mounted on the sample holder and reflectance images were taken continuously every 30 minutes for 24 hours. When studying the effect of muscle stretching, only one side of the muscle was firmly attached to the holder; while the other end was pulled to stretch the whole muscle. Once it was stretched to 125% of its original length, samples were released to the original length. Reflectance images were acquired for each stretched and relaxed position. The whole stretching experiment was finished within 4 hours of the slaughtering.

2.3 Equi-intensity data selection

All acquired images were first smoothed using a 3×3 averaging filter in order to reduce noise. The small area covered by the shadow of the fiber tip was removed from all images. Incident point was considered as the origin. The pixel coordinates of the image were transformed into physical dimensions of the intensity pattern. A pixel corresponding to a specific distance from the origin along the x-axis (along the muscle fiber) was selected. Then all pixels of the same intensity with a ±1% margin were extracted to represent the equi-intensity profile and their coordinates were stored for further processing.

2.4 Evaluation of optical reflectance profile

We designed a non linear fitting model to describe a symmetric geometric shape of equi-intensity contours. A point (x, y) in an equi-intensity contour was represented by:

where a and b describe the axial lengths along the x and y-axis, respectively; θ is the angle with x-axis; and 0 ≤ p ≤ 1. Eqs. (1) and (2) can be simplified to obtain following equation.

or

where $q=\frac{2}{p+1}$ and 1 ≤ q ≤ 2.

Similar to Dagdug et al’s [6] definition of bias parameter for elliptical profile, we defined bias parameter B for our analysis as:

One advantage of this fitting equation is that it can also be used to fit either circles or ellipses. When q =1, Eq. (4) describes a rhombus; while it becomes an ellipse when q=2 or specifically, a circle at q=2 with a=b.

We used the Levenberg-Marquet (LM) nonlinear fitting algorithm [19] to estimate the three fitting parameters a, b, and q in Eq. (4). The fitting errors can be estimated using the average of the algebraic distance from n data points to the fitted curve:

When an image profile was not perfectly aligned with the axes in the laboratory coordinates (Fig. 1), it was rotated to obtain such alignment before applying LM fitting method. In addition, a direct least square error ellipse (DLSEE) fitting method [20] was used initially to calculate the elliptical approximates for each equi-intensity contour. The fitted ellipse center was selected as the origin for the LM method because f(x,y) assumed a coordinate system with origin located in the center of the contour. We multiplied the obtained major and minor axes values of elliptical fitting by 1.2~1.25 and used them as initial guesses in the LM fitting of parameters a and b. It was found that this initial setting provided better stability during convergence. The y value cannot be negative in the fitting because a non integer q value leads to an exception during calculation. Therefore all negative y values were mirrored to the positive side based on axial symmetry. Parameter q was initially set to unity to restrict any overshooting. The a and b values obtained from the LM fitting were then used to calculate the bias parameter B in Eq. (5).

3. Results

Figure 2 shows the optical reflectance and fitting curves for three different samples: an isotropic Intralipid^® phantom, a fibrous extrusion sample, and a Sternomandibularis muscle sample. The muscle sample was acquired from a beef cattle right after the slaughter. Clearly, the reflectance profile in skeletal muscle showed a special pattern which was neither a circle as in Intralipid solution nor an elliptical as in fibrous extrudates.

 

Fig. 2. The reflectance images of (a) 20% Intralipid solution (2.28% concentration); (b) soy extrusion sample; (c) prerigor Sternomandibularis muscle. Images were displayed in pseudo color which represented the pixel intensity in 10-based Logarithm. All images had sizes of 2.54×2.54cm². The graphs in the 2^nd row were the fitting results for equi-intensity contours extracted at 8.0 mm from the incident point. Red dots indicated the extracted raw data points on the contour; while blue lines represented the fitting results. The obtained fitting parameters for the three media were: (a) DLSEE: B=1.01, LM: q=1.98, B=1.02; (b) DLSEE: B=1.40, LM: q=1.94, B=1.42; and (c) DLSEE: B=1.17, LM: q=1.46, B=1.21.

Download Full Size | PDF

Both DLSEE elliptic fitting and Eq. (4) were used to find the best fit for these images. As expected, both methods provided excellent fit for Intralipid^® phantom with bias parameter B close to unity. Fitting using Eq. (4) produced a q=1.98, which indicated a circular profile. Both methods provided an elliptic profile for soy extrudates with B values significantly larger than unity (B=1.40 and B = 1.42). The q value obtained using Eq. (4) was 1.94 (Fig. 2(b)). These results indicated that the light reflectance profile from the fibrous sample was indeed elliptic. However, for skeletal muscle, it is clear that Eq. (4) (blue line in Fig. 2(c)) provided much better fit than DLSEE method (cyan line). As a confirmation of the visual difference, in Fig. 2(c), the error from DLSEE fitting method was 14 times higher than that from the LM fitting method.

 

Fig. 3. The dependence of fitting results on the distance between evaluation location and light incident location.

Download Full Size | PDF

In these experiments, we observed that the shape of light reflectance profiles depended on the distance between the evaluation point and the light incident location. Figure 3(a) shows the obtained B and q parameters at distance from 5 mm to 15 mm for the reflectance image shown in Fig. 2(c). The corresponding fitting errors (Eq. (6)) were plotted in Fig. 3(b). At distances close to the incident point, the fitting results changed rapidly. The calculated bias parameter decreased with the distance. We contributed this phenomenon to the effect of scattering because at larger distance, photons experienced more scattering events that were not related to tissue anisotropic structures. However, starting from ~8 mm along x-axis, the obtained bias parameter changed much slowly. The fitted q parameter was stable between 8 – 13 mm. The calculated fitting error indicated that fitting errors were minimized at distance between 7 – 12 mm. In all subsequent results presented below, we chose 10 mm as the evaluation distance to extract the equi-intensity profiles. However, similar curves were obtained at other evaluation locations.

 

Fig. 4. Change of reflectance profile during stretching. (a) Animated sequence (0.9 MB) of equi-intensity contour (red) with numeric fitting results (blue). The image was 2.47×2.47 cm² in size. [Media 1] (b) Quantitative changes of fitting parameters.

Download Full Size | PDF

Figure 4 shows the effect of stretching on the reflectance profile analyzed using our fitting model (Eq. (4)). A pre-rigor sternomandibularis muscle was stretched to 125% of its original length and then relaxed to its natural length. The measurements were started 48 minutes after the slaughter and repeated four times with 40 minutes interval. The reflectance image acquired for the first relax state was selected as reference. The pixel intensity at 10 mm along the x axis from the origin was used as the reference intensity to extract the contour profiles in all images. Fitting parameters and bias parameter were calculated for each event. Figure 4(a) shows an animated sequence of the changes of the reflectance contour along with the fitted curve during stretch and relax. Figure 4(b) shows the fitted bias parameter B and the q values. It can be seen that the bias parameter was reduced when stretched and recovered when released. Comparing to the B value, the q values had very little change.

 

Fig. 5. Change of reflectance profile during rigor process. (a) Animated image sequence (2.5 MB). The image had a physical size of 3.05×3.05 cm². [Media 2] (b) Fitting parameters q and B versus time. For clarity, data points measured ever hour were shown in the figure.

Download Full Size | PDF

Optical reflectance was also monitored during rigor mortis development at room temperature (Fig. 5). Muscle samples were fixed by restraining both ends to a sample holder. Images were recorded at every 30 minutes. First measurement was taken 35 minutes after the slaughter. Each time, the average intensity of a pixel located 10 mm from the origin was used to extract the equi-intensity contours. Figure 5(a) shows the animated sequence of the reflectance images during the 25 hr rigor process. Figure 5(b) shows the variations of the bias parameter B and the fitting parameter q. These results indicated that the q parameter had little change during the first 8 hrs while the bias parameter B increased from 1.2 to ~1.7. The q then increased until it reached ~1.8 at ~24 hrs after slaughter. On the other hand, the B parameter stabilized from 8 ~ 14 hrs and then started to decrease thereafter. After 24 hrs, the reflectance profiles approached circular patterns.

4. Discussion

Equi-intensity profiles of optical reflectance in Sternomandibularis muscle showed a distinct pattern that can not be explained using currently available anisotropic photon migration models [3–9]. A previous model [4] based on Mie scattering of infinite cylinder can produce similar rhombus patterns. However, in their simulation the pattern became elliptic when the measurement location was only ~2 mm away from the light incident point. On the contrary, the patterns observed in prerigor muscle were distinctly different from ellipse even at locations that were much more than 1 cm away from the incident point (Fig. 2 & 3). We found that such distinct pattern can be quantitatively analyzed using a simple fitting equation (Eq. (4)) with three parameters a, b, and q. As in previously reported studies [6], the non-unity bias parameter B suggested anisotropic optical properties parallel and perpendicular to muscle fibers. More importantly, the smaller q parameter clearly indicated the deviation from the elliptic profile predicted from all existing models.

We found that the significant decrease in bias parameter during stretching was mostly due to the increase in parameter a. This suggested that light propagated longer distances along the muscle fiber. Since stretching pre-rigor muscle resulted in proportional changes in the sarcomere length as well as the fiber diameter [21], both effects had the potential to modulate the photon migration process. On one hand, an increase in sarcomere length can alter the scattering phase function and scatter more light toward fiber orientation [22]. On the other hand, the decrease in fiber diameter may reduce the scattering cross section [15], which reduced the probability that light being scattered away from the fiber orientation. More studies are necessary to elucidate the detailed mechanism.

In rigor mortis analysis, equi-intensity contour were extracted at fixed distance (10 mm) from the incident location in all images. The bias parameter B reflected the relative distance that photon traveled perpendicular and parallel to muscle fibers. Based on our observation in the stretching experiment (Fig. 4), the initial increase in bias parameter during the first 8 hrs seemed to suggest a contracting process of sarcomere. Such sarcomere shortening phenomenon was studied before [23] and found to be subjective to effects of environmental factors especially temperature. As rigor mortis progressed, actin and myosin proteins were bound together and formed actomyosin complex. This stage corresponded the phase in Fig. 4(b) where the bias parameter became stabilized. It has been recognized that after ~12 hrs, the proteolysis process [24] eventually became significant because of the release of autolytic enzymes stored in lysosomes. These enzymes broke down the actomyosin complex and disintegrated sarcomere structure. The destruction of sarcomere structures may explain the decreasing of the bias parameter. More interestingly, the change in q values indicated the shape of equi-intensity profile changed significantly during the rigor process. Initially, q was small and equi-intensity contours formed a pattern closer to a rhombus. After 24 hours, the equi-intensity contour formed a pattern closer to an ellipse as a result of the loss of sarcomere integrity. It is interesting to note that when we kept muscle samples at controlled low temperature (4 °C) during rigor, bias parameter varied slowly with the time due to a slower rigor process. Nevertheless, optical reflectance measurements were able to reflect structural changes of the whole muscle during pre and post rigor.

 

Fig. 6. The reflectance images of (a) prerigor Sternomandibularis muscle acquired in experiments; and (b) Monte Carlo simulation (10⁷ photons). All images were 2.54 × 2.54 cm². The graphs in the 2^nd row show the corresponding equi-intensity contours (red dots) extracted at 10mm from the incident and the LM fitting results (blue lines). Fitting results were (a) q=1.54, B=1.15; and (b) q=1.55, B=1.13.

Download Full Size | PDF

Light diffraction by sarcomere in single fibers is a well studied subject in muscle physiology. As an approximation, muscle fibers are often treated as a volumetric phase grating in theoretical studies. Two dimensional couple wave theory [25] has been successfully applied before in studying muscle sarcomere structures in single muscle fibers or fiber bundles [26]. Our current experimental studies suggested the unique reflectance observed in whole muscle is likely due to the sarcomere diffraction. In order to test this, we applied a three dimensional coupled wave theory [27] in a semi-infinite Monte Carlo model. The sarcomere structure model and properties proposed by Thornhill et al. [28] was used in the calculation. Each sarcomere was taken as a single grating with following geometrical parameters: sarcomere length = 3.1μm, length of A band =1.5μm, and length of I band = 2.0μm. A photon packet experienced regular scattering events as in the conventional Monte Carlo model [28] with Henyey–Greenstein phase function. However, in between two consecutive scatterings, the photon had a probability to be diffracted instead of following a ballistic trajectory. Based on incident angle and the step size, diffraction efficiency and diffraction angle was computed using the coupled wave theory. The effect of step size and incident angle can change the diffraction efficiencies as shown in previous studies [25–28]. The actual diffraction angle was sampled based on the relative efficiencies of all diffraction orders and was used as the new photon direction. In addition, to consider the fibrous effect, different scattering probabilities (scattering coefficients) were assigned to the interaction direction parallel and perpendicular to the muscle fibers [2, 5]. This process was continued until photon lost its energy due to absorption or photon was transmitted out of the medium.

Figure 6(b) shows a simulation result for optical reflectance with a total 30% diffraction probability. A pencil beam was incident upon the thin glass window on top of the sample at 45° oblique incident angle. Optical properties used were scattering coefficient along x-axis μ_sx=13.4cm^-1, scattering coefficient along y-axis μ_sy=30cm^-1, anisotropy g=0.94, absorption coefficient μ_a=0.1cm^-1. The quantitative fitting results are also shown in Fig. 6. It can be seen that the simulated reflectance pattern (Fig. 6(b)) was very similar with experimental results (Fig. 6(b)), especially at larger distance from the incident point. However, discrepancies can be observed at locations that were closer to the light incident location. Figure 7 shows the calculated B and q parameters at different distance for the reflectance images shown in Fig. 6. The simulated images were quite noisy at distances greater that 12 mm and we could not get reliable fitting results. It can be seen that the simulated B and q had very good agreements with experimental observations at distance larger than 8 mm. The differences at smaller distances were obvious. This is likely caused by the aforementioned method we applied to handle the fibrous effects, which may only be valid in the diffuse region. Better results might be produced if the cylindrical geometry [4] of muscle fiber is considered.

 

Fig. 7. The dependence of fitting results (B and q) on the distance between evaluation location and light incident location: comparison between experiment (black) and simulation (red).

Download Full Size | PDF

5. Conclusion

All existing light propagation models failed to describe the light propagation in striated muscles because the effects of sarcomere structure were not accounted for. We have developed a fitting function which can fit the experimental equi-intensity contours very well. The fitting parameters provided a quantitatively description of light reflectance patterns measured at muscle surface. Our in vitro studies on prerigor muscles highlighted that optical reflectance image can monitor the changes in the muscle microstructures. Our Monte Carlo simulation added strong evidence that the sarcomere structure was responsible for the unique reflectance patterns observed in our experiments.

As sarcomere is the fundamental functional unit in muscle [11], a method that can monitor its dynamics in whole muscle may prove useful in many muscle related studies such as exercise physiology. In addition, sarcomere length is an important factor that determines the eating quality of meat [30]. Therefore, such a method may provide a useful tool in characterizing muscle qualities in meat industry.