1. Introduction

Optical propagation of light in biological materials is dominated by scattering due to presence of inhomogeneities in the cellular structures of the order of the incident wavelength [1]. Light propagation within a biological sample can be characterized by 3 optical parameters. These optical properties are denoted by μ_a, μ_s and g. The absorption and the scattering coefficients, μ_a and μ_s, give the probability per unit path length of photon being absorbed or scattered, respectively. The anisotropy coefficient, g, is defined as the average cosine of the photon scattering angle. Different values of the optical parameters give an insight of the structure and component concentration of the biological sample under study. For example, normal human whole blood consist of about 55 vol % plasma (90% water, 10% proteins) and 45 vol % cells (99% red blood cells “erythrocytes”, 1% leukocytes and thrombocyte), giving a high absorption coefficient and a strong forward scattering, due to the big size of the red blood cells [2].

In general, the optical properties of a turbid media are obtained by using different solutions to the radiative transport equation and comparing its predictions with temporal, frequency and/or spatial measurements, being the steady-state intensity measurements the simplest to perform. These solutions usually are approximated and have a limited validity range [3]. For example, Beer’s law neglects scattering and is inappropriate for thick scattering materials. The Kubelka-Munk method and variants are widely used, but are limited in their accuracy [4–9]. Other methods based in the diffusion approximation, like pulsed photothermal radiometry [10], time resolved spectroscopy [11], radial reflectance spectroscopy [12], and an iterative method based in reflection and transmission measurements [13], tend to be more accurate than Kubelka-Munk. However, the diffusion approximation assumes that the internal radiance is nearly isotropic and consequently it is a poor approximation when scattering is comparable to absorption [14]. A more general method, suitable to steady-state intensity measurements, inverse adding doubling (IAD) [15], is an iterative algorithm based in the adding-doubling method [16] and is sufficiently flexible that anisotropic scattering and internal reflection at the boundaries may be included. Unfortunately, IAD is unable to recover the anisotropy coefficient from reflectance and transmittance measurements [15]. In contrast, the Monte Carlo method is flexible and provides an accurate solution to the radiative transport equation based in probabilistic rules for photon propagation [17]. There exists in the literature several Monte Carlo-based inversion algorithms that recover the optical properties of a turbid media from reflectance and transmittance experimental measurements. Some of them use look-up tables of reflectance and transmittance at different optical properties to compare with the experimental data, interpolating to the missing data [18–20], with different integrating sphere’s corrections [19, 21, 22] or no corrections [18, 20], using different methods for the optimization procedure [19, 21, 22], or even hybrid approaches to the light transport model [21, 22].

Optical properties data of blood are important not only for many diagnostic and therapeutic applications in laser medicine but also for routine medical diagnostics. As stated by Meinke et al. optical properties of blood are required for the calculation of the light distribution in blood-perfused tissues, for example, in optical tomography, fluorescence diagnosis, photodynamic therapy, and laser-induced thermotherapy [23]. Due to the importance of knowledge of the optical properties of blood, these properties have been studied from different approaches and using different techniques. The behavior of human blood has been also studied at different hematocrit concentration, temperature, at different incident wavelength, both in oxygenated and deoxygenated blood [2, 17, 23–28].

Human whole blood, in the coagulation process, forms clots due to a chemical change in a soluble protein normally found in the blood, fibrinogen, that makes it insoluble and breaks the ability to interlink with other molecules to form large macromolecular aggregates [29]. It is commonly assume that the formation of clusters can be viewed as a formation of larger particles embedded in the plasma, which can be observed in the behavior of the anisotropy factor, g. For this, a recovery algorithm for g could be sensitive to the formation of clots in human blood.

In this work we present and validate an algorithm to retrieve the optical parameters from steady-state intensity measurements, i.e. total diffuse reflectance, R_d, and diffuse transmittance, T_d, using a system of two integrating spheres. Due to the non-linear effects introduced by the spheres, the experimental values are affected and do not correspond to real reflectance and transmittance measurements. The algorithm presented takes into account this effects by correcting the experimental values measured with the two integrating spheres set-up [30]. The recovery algorithm is based on a combination of the well known Monte Carlo Multi-Layered code, [31, 32] and a genetic algorithm. This new algorithm is called Genetic Algorithms and Monte Carlo Multi-Layered (GA-MCML). We propose MCML as the calculation algorithm of the reflectance and transmittance properties of turbid media from the knowledge of the optical properties of the sample under study. The iterative method used to find the best solution to the inverse problem is based on a genetic algorithm, that simulates the natural evolution of species [33]. We demonstrate that this method is as accurate as IAD, but being able to handle new experimental conditions. In particular, we are interested in the variation of the anisotropy coefficient by changing the anticoagulant concentration (AC) in blood samples. The GA-MCML algorithm for recovery of optical properties is sensitive to small changes in the concentration and size of the cells [34]. The advantages of this approach over other Monte Carlo-based methods are that include geometric information of the sample, allowing us to calculate lost of light, and corrections to the integrating sphere’s measurements, a very efficient search procedure and the lack of the use of look-up tables. Also, GA-MCML is able to recover the anisotropy coefficient from reflectance and transmittance measurements. A disadvantage of GA-MCML method is that this method is computationally intensive.

2. Materials and methods

2.1. Sample preparation

For the GA-MCML validation, several samples were prepared. Four synthetic samples (phantoms) with different optical properties were prepared to take experimental data. The phantoms were made of polyurethane, as a polymeric matrix, zirconium dioxide, as scattering agent and Indian ink as the absorption element. Part A and Part B of polyurethane (WC-781 A/B, BJB Enterprises, Tustin, California, USA) were mixed together with the scattering and absorption elements. De-gassing was carried out in a vacuum pot and the phantom was kept in oven at 60 C for 8 hours. The sample diameter was 6.5 cm with different thicknesses. It was considered that refractive index for all the phantoms are 1.5. Also, we prepare five different colloidal suspensions of polystyrene microspheres (PS02N/6703 and PS03N/5717, Bangs Laboratories Inc., Fishers, Indiana, USA) in pure water as optical parameters standards, at two microspheres diameter (0.39 and 0.59 μm diameter) and different concentrations. The polystyrene microspheres have a small standard deviations in size (around 5%). The samples were put in glass cuvettes of 2 mm thickness. It was considered that the refractive index for these samples is 1.33.

Fresh human blood from a healthy blood donor was placed in eight tubes where different quantities of anticoagulant (Na Citrate) were added with a pipette, the samples were stirred and poured in eight specially designed cuvettes of thicknesses between 0.8 and 1 mm. The measurements were carried out 20 minutes after the extraction of blood at room temperature. The quantity of anticoagulant contained in blood was varied from 0 to 17.5% in volume. The effective refractive index of a tissue is often approximated as the volume-weighted average of the values of its components [35]. Refractive index for whole tissues themselves fall in the range of 1.36 to 1.4 [35]. There are discrepancies between different authors regarding the refractive index of whole blood, thereby we decided to use a volume-weighted average refractive index of plasma (1.33) and erythrocytes (1.48) [23]. Thus the refractive index used for the retrieving of optical parameters was assumed equal to 1.4, which is close to the values reported and used in the literature [26, 28]

2.2. Experimental setup

Experimentally, the diffuse reflectance and transmittance measurements were taken using two identical integrating spheres (model RTS-4Z, Sphere Optics, Uhldingen, Germany) (see [36, 37]). The diameter of each sphere was 100 mm, the entrance port had an aperture of 12.7 mm in diameter, and the sample port had a diameter of 25.4 mm in. The integrating spheres have two silicon detectors attached to the north pole of the integrating sphere and are connected to two identical picoammeters (model 6485, Keithley Instruments Inc, Cleveland, Ohio, USA). The tissue sample was placed between the two integrating spheres. The measurements were performed at 633 nm and at room temperature. The 633 nm wavelength was supplied by a helium-Neon laser head (model HRP120, Thorlabs-Inc, New Jersey, USA). The laser beam diameter at 1/e² of the peak intensity was 0.88 mm, and the beam divergence was about 0.92 mrad. The laser beam was directed into the entrance port of integrating sphere 1, whose exit port was coupled with the entrance port of integrating sphere 2; the sample was situated between the two integrating spheres, as seen in Fig. 1. The exit port of integrating sphere 2 was covered with a cap with a reflective surface identical to the integrating spheres coat.

 

Fig. 1 Experimental setup for R_d and T_d measurements using a system of two integrating spheres.

Download Full Size | PDF

From the current detected in the picoammeters, the total reflectance and transmittance were calculated using the following expressions [38].

2.3. Genetic algorithm and MCML program

GA-MCML is an inverse Monte Carlo program that recovers the optical properties of a turbid media from two experimental measurements: diffuse reflectance and transmittance. GA-MCML compares the simulated diffuse transmittance and reflectance with the experimental ones, and start an iterative process to find the best set of optical parameters that matches the experimental measurements [34, 36, 37]. The optimization technique that uses GA-MCML is a genetic algorithm, and the Monte Carlo code is the well known Monte Carlo Multi-Layered developed by Lihong Wang and Steven Jacques [31, 32].

2.3.1. GA-MCML code

MCML program has been included inside of the GA-MCML code in order to solve the radiative transfer equation through a stochastic procedure. Monte Carlo simulations of photon propagation offer a flexible yet rigorous approach toward photon transport in turbid tissues. The method describes local rules of photon propagation that are expressed as probability distributions that describe the step size of photon movement between sites of photon-tissue interaction, and the angles of deflection in a photon’s trajectory when a scattering event occurs. The simulation can score multiple physical quantities simultaneously. However, the method is statistical in nature and relies on calculating the propagation of a large number of photons by the computer. As a result, this method requires a large amount of computation time [31, 32]. For the purpose of this work, the implementation of MCML as the method for solving the radiative transfer equation is sufficient and straight-forward.

An additional correction is required to consider the light that escapes by the edges of the sample due to its finite size. If this correction is not included, the loss light can be interpreted as absorbed light, giving a higher absorption coefficient [39]. The light loss is a function of sample dimensions, the sample port size in the integrating sphere set-up, the illuminating beam diameter and the optical properties of the sample, and can be as large as 30% [40]. Although difficult to overcome in most of inverse methods, loss of light can be taken into account easily in the MCML program. MCML divides the sample in a finite virtual grid to storage the energy loss in each grid element. Energy loss outside the grid is storage in the closer element of the virtual grid. The size of the grid can be defined by the user. We modified the MCML code to use the virtual grid as a measure of the finite size of the sample and to storage loss light by the edge. If a photon exits the virtual grid, all its energy is stored in a new variable, called lateral escape, and the photon is terminated.

2.3.2. Genetic algorithms

A genetic algorithm is a powerful optimization method based on the mechanism of biological evolution. A genetic algorithm starts with a random set of solutions called population. Each individual in the population, called chromosome, represents a potential solution to the problem. Each chromosome evolves through successive iterations, called generations. In each generation, the chromosomes are evaluated using a fitness function. To create the next generation, an offspring, several operations are applied to all the chromosomes, such as a crossover operator or a mutation operator. The new generation is formed by selecting, according to the fitness values, some of the parents and offspring and rejecting other individuals in order to keep the population size constant [33]. After several generations, the algorithm converges to the best chromosome, i.e. to the individual with the best fitness value. Fig. 2 shows the basic scheme of a genetic algorithm,

 

Fig. 2 The general structure of a genetic algorithm.

Download Full Size | PDF

A genetic algorithm encodes the chromosomes into binary strings. The length of the string depends on the required precision. For example, if the requested domain of variable x_j is [a_j, b_j] and the required precision is q places after the decimal point, the range of the domain of the variable should be divided into at least (b_j − a_j) × 10^q size ranges. The required bits (m_j) for a variable is calculated as follows:

The mapping from a binary string to a real number for variable x_j is computed as follows:

The crossover rate (p_c) is defined as the ratio of the number of offspring produced in each generation to the population size (pop – size). This ratio controls the expected number p_c × pop – size of chromosomes that undergo a crossover operation. Mutation is a background operator which produces spontaneous random changes in some individuals. In a genetic algorithm, mutation serves the role of either replacing the genes lost from the population during the selection process or providing the genes that were not present in the initial population, exploring new possibilities of solutions [33]. The mutation rate (p_m) is defined as the percentage of the total number of genes muted in the whole population and controls the rate at which new genes are introduced into the population for trial. If it is too low, many genes that would have been useful are never tried out; but if it is too high, there will be much random perturbation and the algorithm will lose the ability to learn from the history of the search.

The process of evaluating the fitness of a chromosome consists of convert the chromosome’s encoded solution to its decoded solution, i.e. converting binary string into relative real values ${\mathbf{x}}^k=\left(x_1^k,x_2^k\right)$, k = 1, 2,..., pop – size and to evaluate the fitness function f(x^k). In our case, the fitness function represents the relative error between the experimental measurements and the simulated ones using MCML.

A roulette wheel approach is adopted as the selection procedure of each generation, the best individuals will receive a portion of the wheel greater than the received by the worst, based on the fitness function evaluation. To avoid losing the best solution to the current generation, the program automatically keeps this solution as the first individual of the next generation.

GA-MCML program offers a precision of 3 digits after the decimal point. Each generation has 100 individuals with 28 genes generated randomly. The search range for g was selected to be [0, 1], for μ_s was [0, 200] and for μ_a was [0,10] for the GA-MCML validation and for the blood experiments was [0.95, 1] for g, [0, 4000] for μ_s and [0,100] for μ_a. A Boltzmann function was selected as the fitness function. GA-MCML also includes a cut fitness value. If the best element of a given generation has a fitness value greater than the cut value, the program is terminated. In the case of the validation of the method we use 10,000 photons in the MCML simulation and 30 generations in the genetic algorithm. For the optical properties of blood at different grades of coagulation we use 30,000 photons and 30 generations in the MCML simulation and in the genetic algorithm evolution respectively. In both cases no cut fitness value was used.

2.3.3. Integrating spheres corrections

Using a light propagation model, such as MCML or AD [15], integrating spheres measurements may be used to calculate the optical properties of a turbid media. The effect of the integrating sphere on the measured reflectance has been studied, for example, by Pickering [41, 42]. Unfortunately, inconsistencies in the experimental and theoretical parameters make this method difficult to use [38]. Alternatively, Moffitt developed a mathematical model of integrating sphere measurements for both single sphere and double-sphere arrangements [30]. This derivation is consistent with prior integrating theory work. Briefly, a clear and simple characterization of the effect of the integrating spheres upon measurements is given as a single function for the sphere’s gain. The greatest advantage of the integrating sphere is that it collects nearly all the light reflected or transmitted by a sample, however a significant of this light is collected and redirected to the detector. The exact amount depends on geometric factor as well as the reflectivity of the sphere walls and sample, giving a dependence of the sphere gain on the reflectance of the sample, thus, the sphere gain describe the increase in the irradiance on the detector due to the reflective sphere walls and sample.

The so called, Moffitt’s corrections depends on the geometrical parameters of the integrating spheres used and are included in the latest version of IAD [15, 38]. Also requires the knowledge of the reflectance and transmittance at normal (direct o collimated), $R_d^{\mathit{direct}}$ and $T_d^{\mathit{direct}}$ respectively, and at diffuse incidence, R_d and T_d. The latter would require a complete modification of the MCML code, so, at first approximation, we use $R_d^{\mathit{direct}}\approx R_d$ and $T_d^{\mathit{direct}}\approx T_d$. Moffitt’s corrections were applied to all reflectance and transmittance measurements in this work.

3. GA-MCML validation

As indicated before, we recover the optical parameters of several samples using GA-MCML, and we compare our results with the well-known IAD algorithm [15, 38]. In both cases, Moffitt’s corrections for two integrating spheres were used to correct the experimental reflectance and transmittance. IAD also includes a Monte Carlo subroutine that calculates the light loss by the edge of the sample [38].

The results of the experimental measurements of the synthetic phantoms as well as the solutions of microspheres are shown in Table 1. The first column is the sample ID label, MS039 stands for the solution of microspheres of 0.39μm diameter and MS059 for the microspheres of 0.59μm diameter at different concentrations. PBMS stands for a phantom with scattering agent embedded and PBMSA for a phantom not only with scattering agent embedded but also with an absorption agent. Columns two and three indicate the reflectance and transmittance values experimentally obtained by using a system of two integrating spheres. Column four is the thickness of the sample and the last column shows the anisotropy factor used in the recovery of the optical properties, obtained using Mie theory.

Table 1. Experimental measurements of synthetic phantoms and microspheres solutions used for the GA-MCML validation

View Table | View all tables in this article

From the reflectance and transmittance measurements of Table 1, IAD and GA-MCML were used to calculate the optical parameters of the samples. For the case of the microspheres, Mie calculations were carried out using the concentration of the stock solution. For polyurethane phantoms, Mie calculations are no reliable because of the change of volume due to the polymerization, as well as for the inhomogeneous shape of the scattering agent.

As can be seen from Table 2, IAD and GA-MCML retrieve similar optical parameters, no matter the absorption or scattering coefficient or the thickness and anisotropy factor of the sample. Differences between the scattering coefficient coming from Mie calculations and GA-MCML are between 1.64 and 4.9%. In contrast, the effects of integrating spheres can be quite significant, giving differences between the recovery parameters up to 23% [37]. In the samples with small absorption, a bigger difference between the absorption coefficient retrieved using IAD and GA-MCML was found. The error is attributed to the relatively small amount of photons used in the MCML simulation. However, at a higher absorption coefficient, GA-MCML recovers correctly the optical properties (see sample PBMSA-1 for example).

Table 2. Results of the recovery of the optical properties of synthetic phantoms and microspheres solutions using GA-MCML and IAD. Numbers in parenthesis represents percentage errors relative to Mie calculations

View Table | View all tables in this article

4. Results and discussion

GA-MCML was used to determine the optical parameters of the whole human blood samples at different coagulation grades. The first step was to restrict the solution space for the anisotropy factor (g) to a range from 0.95 to 1 based on previous reported values of this parameter for blood samples [2, 17, 25, 28, 35]. During this step the absorption coefficient μ_a and scattering coefficient μ_s were chosen as the free variables. The goal of this first step was to find a suitable value for μ_a. The second step was to recover the anisotropy factor (g) and the scattering coefficient (μ_s) from two experimentally measured quantities: the diffuse reflectance and the diffuse transmittance using the absorption coefficient (μ_a) recovered from the first step. GA-MCML error was estimated using the mean value and standard deviation of several runs of the GA-MCML program.

The diffuse reflectance and transmittance measurements taken using the two integrating spheres set-up (see Fig. 1), are showed in the Table 3 for eight human blood samples at different anticoagulant concentration, AC.

Table 3. Experimental measurements of diffuse reflectance and diffuse transmittance for whole blood at different anticoagulant concentrations

View Table | View all tables in this article

Table 4 shows the absorption coefficient of human blood at different degrees of clotting, obtained from the procedure explained above. In general, the absorption coefficient tends to decrease by increasing the anticoagulant concentration. It is difficult to compare our results with other reported values because our experiment was performed under different conditions. However, in general, the values of μ_a found in this work for different concentration of anti-coagulant range from 6 to about 16 cm⁻¹ which is in agreement with those reported in the literature [2, 25, 28].

Table 4. Behavior of μ_a of whole blood at different anticoagulant concentrations

View Table | View all tables in this article

The optical parameters of whole blood with at different degrees of coagulation, μ_s and g, were calculated using GA-MCML from diffuse reflectance and transmittance measurements, using the absorption coefficients shown in Table 4. The results are shown in Table 5. This results show that the optical parameters of the human blood are characterized by strong scattering and a strong forward scattering, typical of biological tissues (g>0.7), as found by other authors (see for example [2, 25, 28]).

Table 5. Optical parameters, mu_s and g, of blood at different anticoagulant concentrations

View Table | View all tables in this article

Given the main assumption of this paper about the behavior of blood according to its degree of coagulation, the anisotropy factor g is the key parameter, because it contains information on the clustering of the erythrocytes.

Figure 3 shows the anisotropy coefficient corresponding to the values reported in Table 5 for whole blood. Sodium citrate is commonly used as an anticoagulant for blood. A concentration of 10% is required for a complete lack of coagulation. Around this concentration, the anisotropy factor is around 0.996, value close to those found by Nilsson et al. by computing the anisotropy factor using the T-matrix theory of a red blood cell volume equivalent spheroid [43]. At higher sodium citrate concentration, the anisotropy factor decreases abruptly. The point corresponding to 10% of anticoagulant differs by 0.10% compared with the value for 0% of anticoagulant, whereas this difference is 1.11% respect to the higher sodium citrate concentration. The decrease of the anisotropy factor at higher anticoagulant concentration can be attributed to a change in the refractive index, induced by the presence of the anticoagulant itself.

 

Fig. 3 Anisotropy factor as a function of anticoagulant (sodium citrate) concentration.

Download Full Size | PDF

Due to the formation of clusters, the anisotropy coefficient increases dramatically, which was an expected result [44]. In the area of diagnosis, the fact that blood exhibits high anisotropy factor would indicate an abnormal behavior of red cells and the formation of blood clots.

The scattering coefficient was probably the most difficult parameter to understand, because the presence of anticoagulant depicts an additional volume without scatters, which one might expect would produce a lower value of the scattering coefficient. In addition to this effect, the lack of anticoagulant allows the formation of blood clots resulting in fewer scatters elements, but bigger in size. As shown in Figure 4, μ_s decreases by increasing the concentration of anticoagulant. The points corresponding to 7.5%, 10% and 12.5% had an almost constant scattering coefficient, having an almost linear decreasing at lower and higher anticoagulant concentration (AC). Regarding to the data at 10% of anticoagulant, the first point differs by 38.11% and the final point differs by 213.28%.

 

Fig. 4 Scattering coefficient for human blood as a function of anticoagulant concentration at 633nm wavelength.

Download Full Size | PDF

Human blood has been studied under different experimental conditions, the published results offer a very wide range of values [24, 27, 28, 45]. It is not possible to say whether our results are similar to those of someone in particular. However, the values found of the scattering coefficient are of the order of other published results.

Figure 4 can be divided into two parts: when the blood is clotted (before 10% of anticoagulant) and when the blood is not clotted (after 10%). The second part is consistent with the definition of μ_s, where the scattering coefficient describes a medium containing many scattering particles at a concentration described as a volume density ρ_s. The scattering coefficient is the cross-sectional area per unit volume of medium μ_s = ρ_sσ_s[3]. We are considering σ_s constant, i.e. during the coagulation process, the size of the scattering shadow does not change. But, in the first part, neither σ_s nor ρ_s are constant because the light hits a bigger obstacle by decreasing the anticoagulant concentration (AC) and the total volume is changing. Then, the second part is the only part that is probably linear.

The Figure 5 shows the reduced scattering coefficient for human blood at 633nm calculated from the scattering and anisotropy coefficient. The reduced scattering coefficient was a consequence of the combination of the two quantities plotted in Figures 4 and 5 by relation μ′_s = μ_s(1 − g).

 

Fig. 5 Anticoagulant concentration effect on the reduced scattering coefficient for human blood at 633nm wavelength.

Download Full Size | PDF

5. Conclusions

In this work, a new algorithm for the recovery of optical properties (GA-MCML) of a turbid media is presented and validated. We compare our results with the well know algorithm Inverse Adding Doubling. In general, both methods retrieve similar values, but GA-MCML algorithm is more flexible, because can be adjusted to different geometries and experimental conditions. Only at low absorption coefficients, the GA-MCML method is less accurate, but the precision can be increased using a higher amount of photons used in the MCML simulation. The main disadvantage of the GA-MCML algorithm is the computation time. A complete run of the GA-MCML program takes between 60 and 90 minutes, depending on the optical properties of the sample. The sensitivity of GA-MCML depends mainly in the number of photons used in the MCML simulation. The methodology used in GA-MCML is particularly advantageous for retrieving optical properties of the biological tissues since it does not require an estimated initial value or look-up tables, and does not require complex, thus limited in range, inverse strategies. Also, this technique is sensitive to fine changes in terms in the concentration and size of the scattering elements. GA-MCML contains all the required corrections for a successful recover of optical parameters, as Moffitt’s correction for non-linear effects on reflectance and transmittance measurements using integrating spheres, and information about loss of light by the edges, which can be interpreted as absorption.

Using GA-MCML, we studied the effect of anticoagulant concentration (AC) on the optical properties of human blood. The coagulation effect can be seen as a formation of bigger scattering particles, thus giving a strong forward scattering. Because of this aggregation, the main effect of the coagulation on the optical properties can be seen in the anisotropy factor, having a stronger forward scattering when the blood is coagulated, i.e. when bigger particles are present in the solution. We also found that the retrieved anisotropy factor is consistent with the calculations made using T-Matrix computations. Also, the formation of clusters has an effect on the absorption coefficient, at higher anticoagulant concentration, we found a smaller absorption coefficient.