Introduction

Visible and Near-Infrared light are commonly used in the development of nondestructive diagnostic techniques (e.g. early cancer detection, quality assessment of fruits, vegetables and nuts) [1–3]. The propagation of electromagnetic waves in turbid media is fundamentally described by Maxwell’s equations, which are mathematically rigorous and accurate. However, these have a limited applicability in real biological systems as their solution is considered computationally unfeasible for the microstructural and compositional complexity of biological systems [4,5]. The wave properties of light propagating through a turbid medium tend to be averaged out after several scattering interactions. Consequently, several light propagation models are simplified by ignoring the wave-like behavior. Therefore, modeling light as particles with a quantum energy is considered an acceptable approach for turbid media [4–7].

The radiative transport equation (RTE) provides a mathematical description for the physical phenomenon of energy transfer, in the form of electromagnetic radiation, through turbid media. The energy balance of incoming, outgoing, absorbed and emitted photons of an infinitesimal volume in the medium is taken into account. The medium is represented by the following optical parameters: (1) the absorption coefficient µ_a, (2) the scattering coefficient µ_s and (3) the scattering phase function p(θ). The Henyey-Greenstein phase function is a commonly used approximation for the scattering phase function, represented by a single anisotropy factor g [6–11]. The interaction of light with turbid media can thus be quantified with two important optical phenomena: scattering and absorption.

Since the RTE has no analytical solutions for multiple scattering problems, it must be either approximated or solved numerically [6,7]. Several modeling techniques exist: adding-doubling, Kubelka-Munk theory, four-flux and seven-flux models, discrete ordinate method, diffusion theory and Monte Carlo (MC) methods. Since the latter provide an accurate solution for modeling light propagation in turbid media, they are considered to be the golden standard methodology in light propagation modeling [6]. In MC simulations of light propagation, the movement of photons is discretized and the optical processes are described stochastically using probability density functions [1,2,5–8]. This method assumes macroscopic (bulk) optical properties that extend uniformly over the turbid tissue. As long as the tissue is homogeneous on a macro-scale, this assumption is valid. The expected value of a physical quantity (e.g., reflectance values) is found by running simulations many times over, each time applying the same probability functions, and eventually calculating the average of multiple independent simulations [6]. This makes MC methods perfectly suited for solving problems involving multiple scattering.

Classic MC algorithms, such as the widely accepted MC code for multi-layered media (MCML) [6], are capable of simulating the light propagation in multi-layered, semi-infinite slabs, both in 2D and 3D. This approach assumes that the different plane-parallel layers (e.g., dermis and epidermis of human skin or epidermis and cortex tissue of an apple) are homogeneous media. However, the organization of the different tissue layers can be far from plane parallel. Furthermore, biological tissues have a hierarchical structure consisting of cells with different organelles and intercellular spaces. As their dimensions are in the same order of magnitude as the wavelengths of visible and NIR light, the cell organelles are the dominant scatterers in biological tissues [12].

A typical MC algorithm consists of two distinct parts: physics modeling and geometry tracking. The geometry tracker, also know as the ray-tracer, keeps track of the position and properties of the photon packets that are transported through the medium [13]. A general-purpose code, such as the MCML code, employs simple geometries to track photon movement (planes, cylinders). Some codes have made adjustments by using a voxelized geometry [13,14]. As voxelized images are not well-suited for modeling targets with curved boundaries or locally refined structures, MC methods have been elaborated which can simulate the light propagation in a tetrahedral mesh [13,15]. Prior to applying the MC methodology, this last category of MC methods mesh the voxel based images, derived from tomographic techniques (e.g., X-ray tomography). These mesh-based MC or PenMesh methods provide a good solution for simulating optical profiles through complex structures.

However, in light propagation simulations the physics modeling is at least equally important as the geometry tracking. The scattering phase function is a probability density function which describes the chance that a photon is deflected in function of the scattering angle. The most commonly used example is the Henyey-Greenstein (HG) phase function [1,2,5–8], an analytical function originally derived for modeling scattering by interstellar dust. In the past, the HG phase function has been widely used in biomedical optics because of its convenience and simplicity. Moreover, MC algorithms function computationally more efficient with analytical functions which approximate the shape of the actual phase function, as they require a translation into a cumulative distribution function. As the scattering phase function determines the deflection of the photons at every scattering event, it has a large impact on the simulated results [7]. Therefore, the simplicity of the HG phase function comes at the expense of a reduction in the simulation accuracy [16–20]. In order to promote more accurate simulation of the light propagation through turbid media, the over-simplistic parametric scattering phase functions should be replaced with a more accurate alternative. Turbid media typically consist of a variety of scatterers, each having a specific contribution to the scattering properties of the (uniformous) tissue type in which they are situated. The scattering phase function for different types and sizes of scatterers present in the biological tissue can be calculated with Mie theory (spherical and cilindrical scatterers), T-matrix method (arbitrary geometries), Rayleigh-Gans-Debye approximation (ellipsoids), etc.. The scattering phase function for each uniformous tissue type can be obtained through the combination of the scattering phase functions of the respective scatterers [16–20]. This has already been achieved for a set of plane parallel, homogeneous bulk layers by adapting the commonly used Monte Carlo code for multi-layered tissues [6] to accept any discretized phase function [20]. However, as the assumption of plane parallel, homogeneous bulk layers is not valid for biological tissues, this phase function flexibility should be extended to tissue meshes. Therefore, the main goal of this study was to develop and validate a meshed MC model with flexible phase function choice (MMC-fpf), where the contribution of each type of scatterer can be incorporated. To reach this goal, the Monte Carlo method will be developed for use in combination with the flexible tool to simulate the bulk optical properties of polydisperse suspensions of spherical particles in an absorbing host medium [21]. This combination is essential to allow simulation of the light propagation in microstructured, turbid media such as biological tissues.

The accuracy of the mesh-based MC methods depends on both the functionality of the algorithm and the accuracy of the mesh. A mesh could be generated by simplifying the structure of biological tissues [22]. In that case, however, various assumptions must be made on the shape and size of the cells and on their spatial arrangement in tissues [22]. In particular, cells would be represented through simple geometric objects, and tissues by juxtapositions of such objects [22]. As these approximations often result in a very rough approximation of the real tissue miscrostructure, the aim of this study was to simulate the light propagation through meshes generated from 3D-microstructure images, derived from high-resolution X-ray tomography. The use of such 3D-microstructure images allows to consider the real geometry of the different cell types in the light propagation simulations. However, to be able to perform such simulations, the scattering properties of these cell types should be known. As these are not available in literature and very difficult to measure, the bulk scattering properties of these celly types will be estimated from the size distributions of the different organelles present in these cells. The potential of this approach will be demonstrated for plant leaf tissue, which is commonly modelled as a collection of simplified structures [22–26].

2. Materials and methods

2.1 MMC-fpf algorithm

The MMC-fpf algorithm is a Monte Carlo methodology developed in Microsoft Visual C + + 2010, and shares many properties with other MC codes. Therefore, the description of the algorithm will focus on the differences with a classic Monte Carlo method, the well-known MCML code [6]. The flowchart of the MMC-fpf algorithm is illustrated in Fig. 1. The scheme is similar to that of other Monte Carlo algorithms [6,20], with the important difference of being adjusted to function in a mesh. The photons move in a tetrahedral mesh, defined by 4 nodes and 4 sides [5]. This type of mesh has been chosen in order to approximate curvatures which are present in biological tissues.

 

Fig. 1 Flowchart of the MMC-fpf algorithm.

Download Full Size | PDF

The mechanism of the MMC-fpf algorithm is visualized in 2D in Fig. 2. As can be seen from the figure, the first advantage of a meshed Monte Carlo code is that it can handle curved boundaries by meshing them. The code tracks a photon packet while it propagates through the mesh elements of a specific tissue type and records the distance traveled in each tetrahedron to compute the absortion profile. The scattering is defined by a scattering phase function, which has been computed in advance (e.g. through Mie theory). At the interaction plane between two tetrahedrons with different optical properties, and more specifically different refractive indices, the photon packets can either be refracted or reflected. If the photon packet is refracted into the next tetrahedron, its further propagation is defined by the optical properties which have been assigned to that tetrahedron.

 

Fig. 2 Visualisation of the MMC-fpf algorithm: the figure illustrates how an absorption profile is generated and how the photon packets move from one element to another. The scattering process is defined by a predefined phase function, which can be computed from Mie theory (instead of an Henyey-Greenstein approximation). Even though the resulting anisotropy factor is the same, the phase functions are different.

Download Full Size | PDF

1.1.1. Barycentric coordinates for efficient ray tracing

In order to locate photons as efficiently as possible inside an element, a barycentric coordinate system has been used. This allows to evaluate very fast whether a photon can continue its propagation inside a mesh element. In Fig. 3 the concept of this barycentric coordinate system is illustrated for a triangle (2D):

 

Fig. 3 Barycentric coordinate system.

Download Full Size | PDF

Each point on the triangle can be described with the classical cartesian coordinates. For each element - with nodes [x₁,y₁], [x₂,y₂] and [x₃,y₃] - the coordinate can be rewritten as:

The set of λ-parameters (λ₁,λ₂,λ₃) will form the basis of the barycentric coordinate system. The barycentric coordinates of a point (x,y) are defined as the coordinates of that point relative to the nodes of the triangle with coordinates (x₁,y₁), (x₂,y₂) and (x₃,y₃):

Which is based on the important concept that λ₁ + λ₂ + λ₃ = 1. In matrix notation this becomes:

This matrix formulation can be easily extended to 3 (or more) dimensions. When the position of a traced photon is expressed in barycentric coordiantes, it will be situated inside an element if the following condition is fulfilled:

It should be noted that the T matrix is irrevocably connected to the coordinates of a tetrahedron. This is logical as the coordinates express whether a point (or photon) is situated inside this geometric shape or not. If one of the λ-values is equal to zero, the photon is situated on the boundary of the mesh. The expression of the photon positions in barycentric coordinates thus allows to evaluate very fast whether a photon can continue its propagation inside a mesh element or not.

1.1.2. Movement between elements

At the crossing of the boundary between two different elements, the Fresnel and Snell equations are evaluated to determine whether the photon is transmitted from the first element into the second or reflected back in the first, and to calculate the new directional vector of the photon. The Fresnel and Snell equations are described as a function of the angles between the incident, reflected and refracted rays and the normal of the interface. In a mesh, the boundary is formed by the plane defined by the three mutual nodes of the element and its neighbour. These three points define the equation of the plane, and as a consequence the normal on the plane as well. A rotation matrix is established, which aligns with the z-axis. This rotation matrix can be used to transform the directional vector and thus allows a straightforward computation of the Fresnel and Snell equations. The inverse rotation matrix is then used to calculate the new directional vector, after reflection/refraction.

1.1.3. Defining neighbouring elements

To reduce the computation time of the MMC-fpf algorithm, the 4 neighbouring elements for each element are searched in advance. This allows the algorithm to only evaluate the 4 different neighbours once a photon is no longer situated inside an element, instead of checking all the possible elements. Each neighbor has 3 nodes in common, defining the common plane. If the photon travels a large step with size s, thereby crossing several elements in one step, the pathway is split up at the boundary. The photon first travels to the boundary of it’s current element, ensuring that the boundary interaction is modeled correctly.

1.1.4. Recording absorbance

The MMC-fpf algorithm was developed for simulating light propagation in realistic tissue structures. Therefore, the thetrahedrons of the mesh can potentially be smaller than the mean free path. As a result, a photon could pass through an element without interacting with the tissue. However, a photon passing through that small element, loses a part of its energy due to absorption. Classically, the energy lost due to absorption is only stored at each scattering event [6,20]. The fraction of energy lost ($\Delta W$) at that specific interaction is then attributed to the element or layer where the photon is situated and the photon continues its trajectory with a lower amount of energy.

In order to correctly characterize the absorption, the MMC-fpf algorithm tracks the travelled distance L inside each element. After the interaction at the boundary of the element, the initial weight of a photon W₀ is adjusted to the new weight W using Beer’s law by taking into account the travelled distance L in the medium with absorption coefficient µ_a:

Both formulas are technically interchangeable, with the difference that the first one is a discretized version of the latter. In the discretization process, the travelled distance L is linked to ${\left({\mu }_a+{\mu }_s\right)}^{-1}$.

1.1.5. Mirroring of mesh

The voxel images acquired by 3D microstructure imaging techniques such as X-ray micro-computed tomography (µ-CT) typically have quite limited dimensions because of the fine spatial resolution. Using the meshes obtained from these images directly in the Monte Carlo simulations can be problematic as many photons might reach the boundaries.

Certain thin tissues can be seen as homogeneous along multiple axes. In the case of a plant leaf, the thickness can be captured with the resolution of µ-CT images. Other dimensions, however, are often too extensive to cover in a single mesh with normal computational power. To get around this problem, a mirror is constructed at the edges of the mesh. Photons which are about to exit the tissue through either the top or bottom surface are counted as transmittance and reflectance, respectively. Other planes function as a perfect mirror, making sure the photons propagate inside the volume. The idea is that the meshed structure is a representative subset of the entire tissue. Photons exiting the structure, mirrored back into the tissue, are equivalent to photons entering an identical, neighbouring mesh. This assumption is only valid for light sources that display a radial symmetry with the mesh, as is the case for the uniform light source covering the entire top layer of the mesh, used in the simulations for generating the figures in this paper. In that case, the total simulated reflectance, transmittance and absorbance will be correct. However, if one would be interested to simulate a spatially resolved absorption and/or reflection profile, it should be noted that this mirroring procedure can only be applied for a diffuse light source.

2.2 In silico validation

The meshed Monte Carlo methodology has been validated by comparing the simulation results for a tissue consisting of 4 homogeneous layers with those obtained with MCML [5] (the MCML algorithm models the light propagation through layered tissued, hence the layered meshed tissue), which is considered to be the standard for light propagation modeling in turbid multi layer systems [1,2,5,10]. The goal of this in silico validation is to evaluate the proposed methodology against the widely used MCML code for the same multi-layered sample configurations, while both Monte Carlo algorithms use a different methodology.

The MCML code assumes the tissue to be a slab - a simple multi-layered geometry - and can only take the Henyey-Greenstein phase function into account. Therefore, a tissue mesh has been generated to mimic the multi-layered geometry assumed in the MCML simulations (Fig. 4). A varying range of optical properties has been simulated with both methodologies, and the resulting total reflectance, transmittance and absorbance have been stored for validation purposes. The absorption coefficients were in the 0.01-3.5 cm⁻¹ range to correspond with the typical range of optical properties reported for biological tissues up to 1300 nm [2,3,12,27]. Moreover, the scattering coefficients were varied between 20 and 120 cm⁻¹, while the anisotropy factors ranged from 0.7 to 0.95. Fifty million photons were used in each simulation. The optical configurations used for the in silico validation are summarized in Table 2 in the appendix.

 

Fig. 4 Four-layered homogeneous slab used for in silico validation.

Download Full Size | PDF

The next part of the validation compares the results of simulating the light propagation with a phase function, developed from microstructural information, with its Henyey Greenstein equivalent. In this validation, the microstructure is defined as a collection of spherical particles described by the particle size distribution illustrated in Fig. 5, which is representative for the size distribution of the mitochondria in a leaf tissue. When this particle size distribution is fed to the microscale light propagation tool based on Mie theory described in [21], the phase function illustrated in the bottom left of Fig. 5 is obtained. If we want to approximate this phase function with the Henyey-Greenstein phase function, we should calculate the average cosine of the simulated scattering phase function to define the anisotropy factor g, which in turn defines the shape of the Henyey-Greenstein phase function equivalent. Monte Carlo simulations are performed on the mesh of Fig. 4, with both phase functions to investigate whether an equal average cosine of the scattering phase function is sufficient to obtain the same absorption profile.

 

Fig. 5 (top) Particle size distribution of scatterers in the tissue (bottom left); Phase function computed through Mie theory, starting from the PSD (bottom right); Henyey-Greenstein approximation with the same value of the anisotropy factor.

Download Full Size | PDF

2.3 Experimental validation

2.3.1 Reflectance and transmittance measurements

An experimental validation has been performed on tomato leaves (Solanum lycopersicum L. cv. Growdena). A Double Integrating Sphere (DIS) measurement system was used for measuring the total reflectance (Mr) and total transmittance (Mt) in the wavelength range from 550 to 850 nm, with an interval step of 5 nm [28]. Replicate measurements were carried out on five leaves, each from a different plant of the same cultivar. Only large (older) leaves were used, since a minimal leaf area of 2.5 cm diameter was required for DIS measurements. This set-up has been approximated as a uniform light source in the MMC-fpf simulations. The DIS set-up illuminates the leaf with a collimated beam, which has limited dimensions and is therefore not a uniform light source. However, because of the limited size of the geometrical model of the tomato leaf, the light source was assumed to be uniform over this small mesh.

2.3.2 Geometrical model of tomato leaf

The microstructure of tomato leaves was reconstructed from synchrotron computed laminography experiments that were conducted at beam line ID19 of the European Synchrotron Radiation Facility (ESRF, Grenoble, France) [29,30]. The tomographic reconstruction of a leaflet sample was performed with a filtered back projection algorithm using the PyHST (ESRF) software after correction for sample motion using GNU Octave software [29]. The contrast in the laminography images was insufficient to distinguish cell organelles such as chloroplasts and the vacuole (further details are provided in [29]).

In the computational model, chloroplasts were modeled as a cluster with a thickness of 2.6 µm and a diameter of 6 µm adhering to the mesophyll boundary and occupying about 24% of the modeled mesophyll cell volume [29]. In every cell, a layer of 2.6 µm thickness inside the plasma membrane of the mesophyll cell was first created. The layer was then segmented into individual brick-like patches of similar size (about 6 µm). This procedure provided a voxel-based image representing the complex 3D structure of a tomato leaf, distinguishing between several subcellular structures. The acquired voxel-based 3D microstructure images were converted into 3D tetrahedral finite element (FE) meshes with the ‘Iso2mesh toolbox’ [14]. In such a meshed structure, different optical properties can be assigned to each tetrahedral element in the mesh.

The resulting meshed tomato leaf section is illustrated in Fig. 6. The difference between the palisade mesophyll layer – upper part – and the spongy mesophyl layer – lower part – is clearly visible in this figure. The upper and lower epidermis are visualized as well. The chloroplasts are situated within the mesophyll cells.

 

Fig. 6 (A) Synchrotron computed laminography image of a brick-like section of a tomato leaf; (B) corresponding mesh for MMC-fpf. The different colors are the different types of tissues in the tomato leaf. Different optical properties were attributed to different types of tissue.

Download Full Size | PDF

2.3.3 Optical properties

Knowledge of the optical properties of the different subcellular structures is essential to simulate the light propagation through these tomato leaf meshes. As they cannot be measured directly, they were estimated based on the information available in literature. The absorption properties have been derived from the absorption properties of the main absorbing components reported in literature, whereas the scattering properties have been calculated with Mie scattering theory based on the particle size distribution of the different scattering particles estimated from the values reported in literature [Table 1].

Table 1. Organelles and their sizes used for computing the optical properties of the different subcellular structures

View Table | View all tables in this article

The absorption coefficient of the chloroplasts at 680 nm was calculated from the specific absorption coefficient of chlorophyll of 0.1 cm²/µg [23–26], combined with the average chlorophyll content in tomato leaves of 40 µg cm⁻² [31,32]. The chlorophyll content per cm² of leaf, can be converted to the chlorophyll content per volume unit by correcting for the thickness of the leaf and the volume fraction of the chloroplasts. The absorption coefficients at other wavelengths were calculated relative to the aforementioned values based on the absorption spectrum of chlorophyll reported in [25,26]. The absorption coefficients of the epidermis, vacuole and cytoplasm were set to 1 cm⁻¹ for all wavelengths to account for absorption by cell constituents such as other organelles, proteins and metabolites. This value was calculated by combining the specific absorption coefficient of dry matter (~10⁻³ cm²⋅g⁻¹), that is almost constant in the wavelength region of interest [24–26], with a mass fraction of dry matter of 5%. The air was assumed non-absorbing and the absorption coefficient was set to 0 cm⁻¹.

The scattering of photons in the leaf tissue is dominated by the cell organelles which are of similar size as the photon wavelength. In order to obtain realistic scattering properties, the particle size distribution (PSD) of the most important cell organelles (mitochondria, peroxisomes, nuclei, golgi stacks and ribosomes) in the epidermis and cytoplasm, as well as the grana in the chloroplast layers have been translated into scattering coefficients and phase functions through simulations based on Mie theory [21].

The volume frequency of particles for an organelle i, ($PS{D}_V\left(r\right)$) was assumed to follow a normal probability density function of the particle radius r.

The mean $\overline{r}$ and standard deviation ${\sigma }_r$ were calculated by assuming that the lowest and highest value reported in literature (see Table 1) correspond to the lower and upper 95% confidence limits of a normal distribution [34–38]. For the mesh elements corresponding to the cytoplasm, the contribution of each organelle to the bulk optical properties was computed separately and the final bulk optical properties for the cytoplasm were calculated by summing up these contributions:

where µ_s is the scattering coefficient and $p\left(\theta \right)$ is the phase function. The subscript i indicates the contribution of a significant scatterer, e.g. mitochondria, peroxisomes, nuclei, golgi stack and ribosome-like complexes for the epidermis and cytoplasm, or the grana in the chloroplasts. In Fig. 7 the effect of the different particle size distributions in the different parts in the cell on the resulting phase functions is illustrated, while the effect of the particle size distribution on the scattering coefficients as a function of the wavelength is shown in Fig. 8. As the different subcellular structures in the epidermis cells could not be differentiated in the µ-CT images, it was assumed that the distribution of organelles in the epidermis cells was identical to that in the palisade cells. Therefore, an identical distribution of organelles in a smaller volume was calculated by dividing the original volume fraction by the volume fraction of the cytoplasm.

 

Fig. 7 Resulting phase function at 680 nm, expressed in a logarithmic scale. (A), (B), (C), (D), (E) and (F) correspond to mitochondria, peroxisomes, nuclei, golgi stacks, ribosomes in the epidermis domain, and grana in the chloroplast domain, respectively.

Download Full Size | PDF

 

Fig. 8 Estimated bulk scattering coefficient spectra for the different cell organelles in the epidermis as a function of the wavelength.

Download Full Size | PDF

3. Results and discussion

3.1. In silico validation

In Fig. 9 the reflectance, transmittance and absorbance values simulated with the meshed Monte Carlo methodology for the different optical configurations are plotted against the corresponding values simulated with MCML. One can clearly see the close match between both methodologies with an R² above 0.999. The small differences are most likely caused by the stochastic noise inherent to Monte Carlo simulations.

 

Fig. 9 Scatter plot of reflectance, fraction of absorbed energy and transmittance values simulated with MMC-fpf and MCML for 4 different homogeneous slabs.

Download Full Size | PDF

The absorption profiles simulated with MMC-fpf and MCML for the four-layered tissue are illustrated in Fig. 10. The essential pattern is identical, although some important remarks should be made: (1) A first artefact is the result of translating a mesh based absorption profile, into a voxel based 3D image, and afterwards translating this voxel based 3D image into a 2D image with cylindrical symmetry (illustrated in Fig. 10). Since the tetrahedrons are not uniform in size, the absorption information can be spread out over several voxels, which is especially visible near the source. Furthermore, converting the 3D voxel image into the 2D visualization, results in the occasional artefacts. This is visible as the zones at intermediant radial distances, with larger intensity (Fig. 10, right). These artefacts are not visible in the 3D images, and are introduced as a result of the fact that a larger number of voxels is situated within these radial boundaries. (2) The mirroring process is suited for simulating the total reflectance/absorbance/transmittance, but will result in an overestimation of the spatially resolved absorption profile near the edges of the mesh (in case of a point source). The deeper the absorption profile, the stronger this effect will be. This artefact can be avoided by using a diffuse light source (as in the experimental validation). However, as this is not a standard option in the MCML algorithm, this modification has not been used to simulate the absorption profile illustrated in Fig. 10. In practice, it would, however, be recommended to use this modification when simulating spatially resolved profiles.

 

Fig. 10 Absorption profiles MCML vs. MMC-fpf in log10-scale (fraction absorbed in each pixel, log₁₀ scale).

Download Full Size | PDF

Finally, the particle size distribution presented in Fig. 5 was used to compute the corresponding phase function with Mie theory simulations and its approximation with the Henyey-Greenstein phase function. This resulted in two absorption profiles, similar to the ones presented in Fig. 10. In Fig. 11 the ratio of both absorption profiles (linear scale, no logarithmic conversion) is illustrated. This figure clearly illustrates the effect of the phase function on the simulated absorption profiles. As differences up to a factor 2.5 are observed, it can be concluded that simply using a Henyey-Greenstein phase function with the correct g-value does not provide a valid substitute for simulations with an accurate phase function.

 

Fig. 11 Ratio of absorption profiles generated with MCML vs. MMC-fpf.

Download Full Size | PDF

3.2. Experimental validation on tomato leaves

The total reflectance (M_r) and total transmittance (M_t) spectra acquired for the 5 leaves are illustrated in Fig. 12. A clear absorption peak can be observed at 680 nm as a dip in both the transmittance and reflectance spectra of the tomato leaves. This corresponds to photosystem II, which is also known as P680 [38]. It is clear that the variability between leaves is relatively large (up to 50% for transmittance measurements). This can be attributed to the varying microstructures of the different leaf samples as well as to the varying thickness of the different tomato leaves. It can be observed from Fig. 12 that the simulated reflectance (RMSE = 1.5∙100%) and transmittance values (RMSE = 8.15∙100%) match quite well with the measured ones.

 

Fig. 12 Comparison of measured (solid lines, with confidence intervals) and simulated (stars) transmittance and reflectance spectra for the 5 tomato leaves (Growdena).

Download Full Size | PDF

3.3. Absorption profiles of tomato leaf tissues

In Figs. 13 and 14, the absorption profiles simulated based on the estimated optical properties at 680 nm are illustrated for respectively vertical and horizontal cross-sectional slices of the tomato leaf tissue. In these figures, the absorbance is defined as the fraction of energy stored in each pixel/voxel.

 

Fig. 13 (A) Cross section of a synchrotron computed laminography image (dark blue: air; medium blue: cytosol; light blue: epidermis; red: chloroplast; brown: vacuole); (B) Percentage of the photon energy stored in each pixel, expressed on a log10-scale.

Download Full Size | PDF

 

Fig. 14 (A) Cross section of a synchrotron computed laminography image (dark blue: air; medium blue: cytosol; red: chloroplast; brown: vacuole); (B) Percentage of photon energy stored in each pixel, expressed on a log10-scale.

Download Full Size | PDF

Because of the presence of chlorophyl, the absorption coefficient in the chloroplast is logically higher than that of the surrounding cell organelles. Therefore, the absorbed fraction of light is much higher, ensuring an efficient photosynthesis. Furhermore, it can be seen in Fig. 13 that the fraction absorbed in the upper part of the leaf is higher compared to the lower part, due to the location of the light source. The difference is especially noticeable when comparing the upper and lower epidermis layers.

In the case of the horizontal cross-section (Fig. 14) no clear effect of the position on the cross-section can be observed, because the energy is provided by a diffuse light source, emitting photons uniformally over the top layer. However, the absorption differences in the chloroplasts are again clearly visible.

Another way to present the data is by making a 3D figure of the different tissues. It is difficult to control the validity of the 3D absorption profiles quantitatively as no measurement technique is available for this. Therefore, a qualitative evaluation of the meshed Monte Carlo simulations can be achieved by interpreting the absorption profiles.

In Fig. 15 the simulated absorption profile of the cytosol layers at 1200 nm is illustrated, to demonstrate the potential in the NIR spectrum. Absorption in this layer is strongly influenced by the presence of water, which is a strong absorber at this wavelength. The light source in this simulation is a point source, emitting photons from the center of the top layer. The goal is to illustrate the possibility of simulating the light propagation starting from a point source, wich is only possible since the simulation results no longer need to be validated with measurements. The absorption profile has a different radial symmetry compared to Fig. 13 and Fig. 14. This adjustment in the definition of the light source was chosen to demonstrate the flexibility of the MMC-fpf algorithm. The absorbed fraction of energy increases as the cytosol is more closely situated to the light source. This 3D figure provides an immediate overview of the absorption profile in this tomato leaf extract, but provides less information on the absorption at different depths than Fig. 13, as the majority of the cell organelles are hidden behind the outer structures.

 

Fig. 15 Simulated 3D absorbance profile (percentage of energy absorbed) for the cytosol layers where absorbance is expressed on a log10 scale.

Download Full Size | PDF

3.4. Discussion

The MMC-fpf has been succesfully validated experimentally, starting from the microstructural information provided by synchrotron X-ray computed laminography. This approach is a significant improvement in flexibility from modeling light propagation in simplified biological tissues – using either a multi-layered approach or simplistic geometric objects.

Although both the in silico and experimental validation were successful, it should be noted that the approach followed in this study has some limitations which should be tackled in future research. First, the µ-CT images of the tomato leaves were limited in size, such that they only cover the microstructure of a small volume of the leaf which had to be mirrored to obtain a sufficiently large mesh for the light propagation simulations. However, this approach is only valid for diffuse light sources. As a consequence, the pencil beam of the DIS set-up ( ± 2-5 mm diameter) had to be approximated as a diffuse light source. The simulations of light propagation in tomato leaves, demonstrated in this paper, are therefore only valid for normal illumination, because in vivo tomato leaves do not only obtain light directly from one single light source (e.g., the sun). Photons can also be reflected/scattered/transmitted by different neighbouring objects, changing the lighting pattern of a single tomato leaf. So, this complex illumination should be taken into account when simulating the absorption profile in practical conditions. As the conditions for the mirroring are no longer fulfilled in that case, the microstructure of a larger volume should be measured and meshed for such simulations. However, as this voxel-based image has to be converted into a mesh, an overly complex/large tissue might not be fully covered with the meshing procedure.

Another point of attention are the limited spatial dimensions of the different tissue types in leaves. Monte Carlo simulations are a useful tool for geometrical optics, describing light propagation in terms of rays. The photons propagate in a rectilinear path as they travel in a homogeneous medium. Deflections from this straightforward movement can occur at the interfaces of dissimilar media. This is a decent approximation when the wavelength is small compared to the size of the medium it interacts with.

The accuracy of this methodology can be improved by modifying the phase function [20,21]. Therefore, future research is needed to extend this code and include other types of scatterers (non-spherical) present in biological tissues. An accurate phase function makes it possible to take the effect of small particles into account. This is especially necessary when one is interested in spatially resolved optical information close to the light source (< 10 mean free paths). This approach is valid, as long as the scattering medium is homogeneous in each tissue type. In comparison to most other MC codes, the MMC-fpf algorithm is not restricted to (infinite) multi-layered tissues. Thanks to the meshing procedure, it becomes possible to freely define the homogeneous structures in which the photons propagate. However, because of the limited size of these homogeneous structures in tomato leaves, the mesh operates on a small scale and as a result, individual scattering particles (cells, chloroplasts) may be incorrectly modeled as uniformous structures. In a uniform structure, the scattering is caused by an average particle, which is the result of the contribution of all types of scatterers present in the tissue. This approach is unsuited for small structures (the magnitude of a chloroplast is comparable to one mean free path), since the interaction in these structures is no longer averaged out by the large number of scattering events. At this scale, one can theoretically not ignore the wave characteristics of light. However, the Monte Carlo methodology operates by simplifying reality and ignoring the wave-like behavior. Ideally, one would take the wave properties of light into account by describing the propagation through solving the Maxwell equations. This would, however, require other assumptions and lead to a very high computational complexity. While the experimental validation in this study was successful, it is important to keep this in mind when meshing small, complex structures. The experimental validation in this study only involved the total reflectance and transmittance, which are less sensitive to the correct modeling of scattering inside the leaf tissue. However, the absorption profile will be more sensitive to this problem. In theory, the spatially resolved reflectance or transmittance profile could be measured, and compared to the simulation results. However, the spatial scale of such a measurement set-up is typically much larger than the dimensions of the mesh, which might result in a mismatch between the measurements and the simulations.

From the derived bulk optical properties, it can be deducted that the mean free path (the average distance traveled by a photon between successive interactions with the tissue) in the cells is in the range of 30-50 µm. This means that the methodology only provides useful information after the simulated photons have propagated in a volume with dimensions in the range of 30 – 50 µm: (1) The cell and vacuole layer are in agreement with this restriction, while only the thickness of the cytosol and the chloroplast layer of the cells are too low to satisfy this requirement. Furthermore, (2) the limit of one mean free path is a rule of thumb, but not a black and white border of validity. The main motivation for working with a tool for deriving the optical properties obtained from microstructural information was to extend its range of validity. Since it is difficult to quantify the impact of the limited thickness on the accuracy of the methodology, we made a note on it and suggest future research to focus on this quantification. (3) The simulation results gave a good match with the experimental validation. So, while these limitations may have led to some inaccuracy, this did not lead to erroneous results. Finally, (4) it should be noted that the alternative for using the radiative transfer theory would be to solve Maxwell’s equations over the mesh. Although this might be appealing from a theoretical point of view, this would require many other assumptions and lead to a very high computational complexity which makes that this is still far from feasible in practice for complex systems such as plant tissue. Even though this methodology, as any other methodology, has its limitations, it is a substantial step forward compared to the methods typically used to simulate light propagation in plant tissue.

The computation times (performed on a AMD Phantom II X6 1100T Processor, 3.30 GHz) for the simulations presented in this manuscript provide a good basis for estimating the difference in computation time between the proposed algorithm and the MCML code. Roughly speaking, the simulation times increase with a factor if10-100. For example, the results of the in silico validation were obtained by simulating the path of 50 milion photons, both with MCML as with MMC-fpf. The computation times required on average respectively 0.36 hours and 14.8 hours (increase of 42.9). However, the maximum mesh element size was restricted, to obtain an accurate resolution of the absorption profile. Not requiring this fine mesh resolution, reduced the total computation time to 4.231 hours, which still presents an increase in computation time of 12.21 compared to the MCML algorithm. The propagation methodology is more complex, because the mesh boundary interactions have to be taken into account. As a consequence, the computation time is largely influenced by the resolution of the mesh. The smaller the mesh elements, the more computation power is needed to simulate the propagation of light. Another important aspect for comparing the computation times of both methodologies, is the absence of cylindrical symmetry in the MMC-fpf algorithm. As a consequence, a larger number of photons (~10 times more) are necessary to converge to simulations with a similar noise-profile.

4. Conclusions

A meshed Monte Carlo code with free phase function choice has been elaborated. It has a similar approach as applied by other Monte Carlo methodologies with the main differences that it allows a free phase function choice and is capable of tracing photon trajectories through a meshed structure of arbitrary complexity – opposed to working with homogeneous bulk layers. This code has been successfully validated in silico by comparing the simulated quantities for a meshed collection of homogeneous bulk layers to those obtained with the MCML code. The coefficients of determination between the different reflectance values, absorbance values and transmittance values were all above 0.999, which indicates a very good match. The MMC-fpf algorithm has also been experimentally validated by comparing the measured and simulated Vis/NIR total reflectance and transmittance values for tomato leaves. The microstructure of these tomato leaves has been measured by means of synchrotron X-ray computed laminography, allowing an actual 3D microstructure to be used in a light propagation model. Considering that the optical properties of the different tissue components had to be estimated, the RSME values, for both total reflectance (RMSE = 1.5∙100%) and transmittance (RMSE = 8.15∙100%) were very good.

The proposed meshed Monte Carlo method with free phase function has thus proven to be a good and accurate solution for simulating light propagation through biological tissues characterized by a complex 3D microstructure consisting of multiple components with different optical properties. Therefore, this code allows one to perform more realistic simulations of the light propagation in microstructured materials (like tissues) compared to the commonly used codes which are restricted to homogeneous bulk layers.

5 Appendix

Table 2.

Table 2. Overview of the optical configurations used for the in silico validation. All the simulations use the same meshed multi-layered grid. Simulations for 1-layered tissues were accomplished by using the same optical properties for all 4 layers of the mesh.

View Table | View all tables in this article

Acknowledgments

The authors gratefully acknowledge the Institute for the Promotion of Innovation through Science and Technology in Flanders (IWT-Flanders) for the financial support through the GlucoSens project (SB-090053) and the PhD grants of Rodrigo Watté and Robbe Van Beers (grant numbers 101552 and 131777). Ben Aernouts has been funded as a Ph. D. fellow of the Research Foundation-Flanders (FWO, grant 11A4813N).

References and links

1. V. V. Tuchin, Tissue Optics: Light Scattering Methods and Instruments for Medical Diagnosis, 2nd ed. (SPIE, 2007), 840.

2. V. V. Tuchin, Handbook of Optical Sensing of Glucose in Biological Fluids and Tissues (Chemical Rubber Company, 2008).

3. B. Nicolaï, K. Beullens, E. Bobelyn, A. Peirs, W. Saeys, K. Theron, and J. Lammertyn, “Nondestructive measurement of fruit and vegetable quality by means of NIR spectroscopy: A review,” Postharvest Biol. Technol. 46(2), 99–118 (2007). [CrossRef]  

4. M. Kumar, Quantum: Einstein, Bohr, and the Great Debate about the Nature of Reality (W. W. Norton & Company, 2011).

5. A. Ishimaru, Wave Propagation and Scattering in Random Media, Volume 1: Single Scattering and Transport Theory, 1st ed. (Academic, 1978), 250.

6. L. Wang, S. L. Jacques, and L. Zheng, “MCML--Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [CrossRef]   [PubMed]  

7. A. Kienle and M. S. Patterson, “Determination of the optical properties of turbid media from a single Monte Carlo simulation,” Phys. Med. Biol. 41(10), 2221–2227 (1996). [CrossRef]   [PubMed]  

8. A. Welch and M. Van Gemert, Optical-Thermal Response of Laser Irradiated Tissue (Plenum, 1995).

9. T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. 19(4), 879–888 (1992). [CrossRef]   [PubMed]  

10. F. Bevilacqua and C. Depeursinge, “Monte Carlo study of diffuse reflectance at source–detector separations close to one transport mean free path,” JOSA A 16(12), 2935–2945 (1999). [CrossRef]  

11. C. K. Hayakawa, J. Spanier, F. Bevilacqua, A. K. Dunn, J. S. You, B. J. Tromberg, and V. Venugopalan, “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt. Lett. 26(17), 1335–1337 (2001). [CrossRef]   [PubMed]  

12. J. R. Mourant, T. Fuselier, J. Boyer, T. M. Johnson, and I. J. Bigio, “Predictions and measurements of scattering and absorption over broad wavelength ranges in tissue phantoms,” Appl. Opt. 36(4), 949–957 (1997). [CrossRef]   [PubMed]  

13. A. Badal, I. Kyprianou, D. P. Banh, A. Badano, and J. Sempau, “penMesh--Monte Carlo radiation transport simulation in a triangle mesh geometry,” IEEE Trans. Med. Imaging 28(12), 1894–1901 (2009). [CrossRef]   [PubMed]  

14. A. Badal and A. Badano, “Accelerating Monte Carlo simulations of photon transport in a voxelized geometry using a massively parallel graphics processing unit,” Med. Phys. 36(11), 4878–4880 (2009). [CrossRef]   [PubMed]  

15. Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express 1(1), 165–175 (2010). [CrossRef]   [PubMed]  

16. G. Bernard, E. Tinet, J. M. Tualle, and S. Avrillier, “Phase function simulation in tissue phantoms : a fractal approach Phase function simulation in tissue phantoms : a fractal approach,” Pure Appl. Opt. 377 (1996).

17. S. K. Sharma and S. Banerjee, “Volume concentration and size dependence of diffuse reflectance in a fractal soft tissue model,” Med. Phys. 32(6), 1767–1774 (2005). [CrossRef]   [PubMed]  

18. S. K. Sharma, S. Banerjee, and M. K. Yadav, “Light propagation in a fractal tissue model: a critical study of the phase function,” J. Opt. A, Pure Appl. Opt. 9(1), 49–55 (2007). [CrossRef]  

19. D. Passos, J. C. Hebden, P. N. Pinto, and R. Guerra, “Tissue phantom for optical diagnostics based on a suspension of microspheres with a fractal size distribution,” J. Biomed. Opt. 10(6), 064036 (2005). [CrossRef]   [PubMed]  

20. R. Watté, B. Aernouts, and W. Saeys, “A multilayer Monte Carlo method with free phase function choice,” Proc. SPIE 8429, 84290S (2012).

21. B. Aernouts, R. Watté, R. Van Beers, F. Delport, M. Merchiers, J. De Block, J. Lammertyn, and W. Saeys, “A flexible tool for simulating the bulk optical properties of polydisperse suspensions of spherical particles in an absorbing host medium: experimental validation,” Opt. Express 22, 20112–20238 (2014). [CrossRef]  

22. G. M. Palmer, N. Ramanujam, M. M. Verstraete, and S. L. Ustin, “Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Theory and validation on synthetic phantoms,” Appl. Opt. 45(5), 1062–1071 (2006). [CrossRef]   [PubMed]  

23. J. B. Feret, C. François, G. P. Asner, A. A. Gitelson, R. E. Martin, L. P. R. Bidel, S. L. Ustin, G. le Maire, and S. Jacquemoud, “PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments,” Remote Sens. Environ. 112(6), 3030–3043 (2008). [CrossRef]  

24. J. B. Féret, C. François, A. Gitelson, G. P. Asner, K. M. Barry, C. Panigada, A. D. Richardson, and S. Jacquemoud, “Optimizing spectral indices and chemometric analysis of leaf chemical properties using radiative transfer modeling,” Remote Sens. Environ. 115(10), 2742–2750 (2011). [CrossRef]  

25. J. B. Féret, C. François, G. P. Asner, A. A. Gitelson, R. E. Martin, L. P. R. Bidel, S. L. Ustin, G. le Maire, and S. Jacquemoud, “PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments,” Remote Sens. Environ. 112(6), 3030–3043 (2008). [CrossRef]  

26. S. Jacquemoud, “Comparison of four radiative transfer models to simulate plant canopies reflectance direct and inverse mode,” Remote Sens. Environ. 74(3), 471–481 (2000). [CrossRef]  

27. A. N. Bashkatov, E. A. Genina, V. I. Kochubey, and V. V. Tuchin, “Optical properties of human skin, subcutaneous and mucous tissues in the wavelength range from 400 to 2000 nm,” J. Phys. D Appl. Phys. 38(15), 2543–2555 (2005). [CrossRef]  

28. B. Aernouts, E. Zamora-Rojas, R. Van Beers, R. Watté, L. Wang, M. Tsuta, J. Lammertyn, and W. Saeys, “Supercontinuum laser based optical characterization of Intralipid® phantoms in the 500-2250 nm range,” Opt. Express 21(26), 32450–32467 (2013). [CrossRef]   [PubMed]  

29. P. Verboven, E. Herremans, L. Helfen, Q. T. Ho, M. Abera, T. Baumbach, M. Wevers, and B. M. Nicolaï, “Synchrotron X-ray computed laminography of the 3-D anatomy of tomato leaves,” Plant J. (to be published).

30. E. Herremans, P. Verboven, T. Defraeye, S. Rogge, Q. Ho, M. Hertog, B. Verlinden, E. Bongaers, M. Wevers, and B. Nicolai, “X-ray CT for quantitative food microstructure engineering: The apple case,” Nucl. Instr. & Meth. Phys. B 324, 88–94 (2014). [CrossRef]  

31. C. Daughtry, “Estimating corn leaf chlorophyll concentration from leaf and canopy reflectance,” Remote Sens. Environ. 74(2), 229–239 (2000). [CrossRef]  

32. M. Steele, A. A. Gitelson, and D. Rundquist, “Nondestructive estimation of leaf chlorophyll content in grapes,” Am. J. Enol. Viticult 59, 299–305 (2008).

33. G. P. Henderson, L. Gan, and G. J. Jensen, “3-D ultrastructure of O. tauri: electron cryotomography of an entire eukaryotic cell,” PLoS ONE 2(8), e749 (2007). [CrossRef]   [PubMed]  

34. J. R. Austin 2nd and L. A. Staehelin, “Three-dimensional architecture of grana and stroma thylakoids of higher plants as determined by electron tomography,” Plant Physiol. 155(4), 1601–1611 (2011). [CrossRef]   [PubMed]  

35. W. Hong, “SNAREs and traffic,” Biochim. Biophys. Acta 1744(2), 120–144 (2005). [CrossRef]   [PubMed]  

36. A. Verschoor, J. R. Warner, S. Srivastava, R. A. Grassucci, and J. Frank, “Three-dimensional structure of the yeast ribosome,” Nucleic Acids Res. 26(2), 655–661 (1998). [CrossRef]   [PubMed]  

37. T. A. Dittmer, N. J. Stacey, K. Sugimoto-Shirasu, and E. J. Richards, “LITTLE NUCLEI genes affecting nuclear morphology in Arabidopsis thaliana,” Plant Cell 19(9), 2793–2803 (2007). [CrossRef]   [PubMed]  

38. G. Raszewski, B. A. Diner, E. Schlodder, and T. Renger, “Spectroscopic properties of reaction center pigments in photosystem II core complexes: revision of the multimer model,” Biophys. J. 95(1), 105–119 (2008). [CrossRef]   [PubMed]