Additive Approach to Simulation of Malignant Neoplasms Using the Monte Carlo Method

The paper is devoted to additive simulation of Raman light scattering by skin cancer using the Monte Carlo method. Raman light scattering from normal skin, malignant melanoma and basal cell carcinoma is investigated. Based on the photon transport algorithm proposed by L. Wang and S. L. Jacques, a two-stage algorithm for simulating Raman light scattering from skin has been developed. A method for additive modeling of skin pathologies is proposed. The main idea of this method is a hypothesis that an experimental Raman spectrum of normal skin, obtained by averaging in vivo Raman spectra of normal skin, may be served as a “substrate” for the feature simulated Raman spectrum. Thus, the pathology, for their part, may be “grown” by adding on this “substrate” Raman specific components set related to a tumor type. Additive simulation of malignant melanoma on various stages and basal cell carcinoma has been carried out. The possibility of using the developed algorithm to determine the component composition of the skin by the in vivo Raman spectrum of skin is discussed. An attempt to evaluate the change in the concentration of skin components during the development of cancer has been made. © 2020 Journal of Biomedical Photonics & Engineering.


Introduction
The variety and prevalence of various types of skin diseases in the Russian Federation makes us look for diagnostic methods that can provide high accuracy in determining the specific type of newly formed tissue and provide information about its size and topology. In recent years, optical methods such as fluorescence analysis [1], backscattering spectroscopy [2], and Raman spectroscopy [3] have increasingly been used to diagnose and study malignant neoplasms.
Raman scattering spectroscopy allows one to obtain information on the features of the local structure of skin neoplasms by spectral samples of substances in the biological tissue. Since these spectral samples are a set of narrow peaks, the identification of which in the general spectrum is not difficult, the Raman spectra of various tissues are specific and can be used to successfully differentiate various pathologies in tissues [3]. Due to its optical properties, skin tissue is ideal for applying this optical method.
Recent technical advances in the development of equipment for Raman spectroscopy and probe designs allow their implementation in clinical applications [4,5]. The possibility of in vivo measurements, noninvasiveness and accuracy make optical methods, in particular Raman spectroscopy, the undoubted leaders among modern methods of analysis of skin lesions [6].
The registered Raman signal is formed as a result of scattering of light by the heterogeneous structure and J of Biomedical Photonics & Eng 6 (3) 30 Sep 2020 © J-BPE 030302-2 biochemical components of the skin tissue [7]. A comparative analysis of the Raman spectra of malignant and benign neoplasms allows us to determine specific types of skin tissue. In the process of development of a pathological neoplasm, changes in the biochemical composition of skin tissue are occurred that leads to a change in the Raman spectra [6,8]. This makes it possible to determine the stage of the disease and evaluate changes in the composition of the skin by changes in the registered spectra.
Therefore, it is important to have an accurate model to illustrate Raman scattering of light in normal skin and skin neoplasms.

Skin model
We have used a two-layer skin model contained epidermis and dermis (see Fig. 1, top left). To execute a software algorithm, it is necessary to have a discrete model of the tissue. Therefore, the skin model is presented as a cube of voxels. Each voxel is a particular type of tissue and it is described by parameters such as absorption coefficient, scattering coefficient, anisotropy factor, and refractive index [9,10].
The size of the model is 3×3×5 mm, i.e. 300×300×500 voxels. The epidermis layer is modeled equal to 0.2 mm, and the dermis layer is 4.8 mm. The main difference between epidermis and dermis in our model lies in different optical parameters mentioned above (see Table 1). They have an impact on the absorption of light in the tissue volume, which in turn affects the level of the Raman scattering. The light receiver is the upper surface of the skin, the registration of scattered light is simulated. In other words, getting on the plane z = 0 from the tissue volume, the photon is considered captured by the detector.
We have used the approach described in Ref. [9] to determine the optical characteristics of the skin model. Some data for the calculation are taken from Ref. [10]. An example of the calculated optical characteristics for the proposed skin model is presented in Table 1. In addition to these "background" optical parameters, each voxel is characterized by the content of skin components that form heterogeneities in the tissue due to pathological processes. A number of skin basic components were identified in Ref. [7], the change in the content of which characterizes some pathological changes in the skin: collagen, elastin, keratin, cell nuclei, triolein, ceramides, melanin, and water. These Raman active components were used in the developed optical model of the skin.
When modeling basal cell carcinoma (BCC), no changes were made to the geometric structure of the optical model. When modeling malignant melanoma (MM), the tumor was modeled as a cylinder with a diameter of 2 mm, located perpendicular to z = 0, with a melanin content of 30% (see Fig. 1, top right). We used Breslow classification of MM stages: the stage (i.e. tumor thickness) can be changed by changing the height of the cylinder.
Thus, the type of layer (and "background" optical parameters) and the content of the Raman active components in the layer are set for each layer of the optical model. The number of layers of the optical model can be increased. Moreover, if it is necessary to simulate a complex geometric structure, one can set custom optical properties for each model voxel.

Additive simulation approach
The main idea of the method of additive simulation is a hypothesis that an experimental Raman spectrum of normal skin, obtained by averaging in vivo Raman spectra of normal skin, may be served as a "substrate" for the feature simulated Raman spectrum. Thus, the pathology, for their part, may be "grown" by adding on this "substrate" Raman specific components set related to a tumor type. Fig. 1 shows an illustration of the proposed approach. Melanin has been chosen as the main component of the skin, which forms heterogeneities in the tissue caused by pathological processes in MM. Thus, averaged in vivo Raman spectra of normal skin serves as a substrate, and the MM may be grown by adding on this substrate melanin.

Clinical data description
We use the averaged Raman spectra of normal skin and skin with MM and BCC. These spectra were obtained by averaging in vivo Raman spectra of the skin of several patients. In vivo Raman spectra were registered using a portable setup (light source with 785 ± 0.1 nm central wavelength) [6].
In the present study, 34 spectra of skin neoplasm have been used: 24 MM (II stage -9, III stage -13, IV stage -2) and 10 BCC. For each patient, Raman spectra of normal skin and skin neoplasm were registered. The registered spectra of neoplasms were preprocessed with baseline removal, smoothing by the Savitzky-Golay method, data normalization and centering [6]. Some averaged Raman spectra used in this work, are presented in Fig. 2.
All in vivo studies were conducted on patients older than 18 years and with their consent. The studies were

Monte Carlo algorithm
The Monte Carlo (MC) method is a group of methods for solving mathematical problems by modeling random variables. From the point of view of solving the radiation transfer equation, the MC method consists in computer simulation of random trajectories of N of photons [11,12].
Modeling of Raman scattering using the MC method is described many times, for example in Refs. [13,14,15,16]. Matousek et al. [17] have proposed the first elementary model to illustrate the variation of Raman scattering in tissue volume using the MC method. Reble et al. [18] using MC simulation of Raman scattering in a two-layer skin model, have calculated depth-dependent sensitivity of probes based on a single fiber and multifiber. Mo et al. [19] have proposed a method for simulating Raman scattering of light by the MC method for an endoscopic fiber-optic probe. Wang et al. [20] have considered the simulation of Raman scattering by the MC method in multilayer tissue and obtained Raman scattering spectra of normal skin tissue. Nevertheless, the mentioned works do not consider Raman scattering of light in biological tissue with significant heterogeneities caused by pathological processes (skin neoplasms).
The proposed algorithm consists of the main program (see Fig. 3) and the photon transport subprogram developed on the basis of the algorithm proposed by L. Wang and S. L. Jacques [21].
In the first step, the program simulates the propagation of the incident photons through the sample, which results in a distribution of the excitation photons program is executed sequentially at each wavelength from some interval, for which it is necessary to obtain a Raman spectrum. In other words, the same subroutine by L. Wang and S. L. Jacques is executed, but only one of the frequencies from the above interval is taken as the established frequency, and a point source is modeled as the isotropic source. A feature of the developed algorithm is a method for setting the initial weight of the photon packet W R , which allows taking into account the spatial structure of Raman scattering: where F ex (x, y, z) is an element of a three-dimensional array of absorbed radiation which represents results in a distribution of the excitation photons at the first step of the program, µ R (λ) is comparable with the Raman probability and depends on the frequency of the emission photons (see Eq. 2), S R1 (λ) +…+ S RN (λ) are the intensities of the Raman spectra of pure 1…N skin components at the wavelength λ, C 1 (x, y, z)…C N (x, y, z) are the concentrations of 1…N components in the voxel with coordinates x, y, z.
In Eq. (1), µ R (λ) is determined as in Ref. [22]: As a result, a series of two-dimensional arrays are obtained by the number of wavelengths from the above interval. For each array, a summation of all its elements is carried out. The obtained value is the value of the Raman spectrum of the model at the corresponding wavelength. The normalized sums represent the Raman spectrum.

MC simulation of MM of various stages
In this experiment, the light source has been modeled as a collimated beam with a frequency of 785 nm located in the middle of the z = 0. The direction of light is perpendicular to the volume of the tissue. Part of the Raman spectrum of the skin has been simulated in the range from 1200 to 1800 1/cm at 10 1/cm intervals. The geometric structure of the optical model is shown in Fig. 1, top right. Melanin has been chosen as the main component of the skin, which forms heterogeneities in the tissue caused by pathological processes in MM [8]. So, averaged in vivo Raman spectra of normal skin serves as a substrate, and the MM may be grown by adding on this substrate melanin. As a result, Raman spectra of the skin with MM of various thicknesses were obtained (see Fig. 4b). As can be seen from Fig. 4b, with an increase in the tumor thickness h, the Raman spectrum of skin changes predominantly in the ranges from 1325 to 1425 1/cm and from 1500 to 1650 1/cm where the melanin spectrum shows large intensity values in the indicated ranges.
The Raman spectra obtained using the algorithm have been compared with the averaged in vivo Raman spectra of normal skin and skin with various stages of MM. In the present study, 24 spectra of MM were registered using the portable spectroscopic system: II stage -9, III stage -13, IV stage -2. For each patient, Raman spectra of normal skin and skin with MM were recorded. The averaged Raman spectra are presented in Fig. 4a.
Both in vivo spectra and spectra obtained using the algorithm demonstrate an increase in intensity values in the ranges from 1325 to 1425 1/cm and from 1500 to 1650 1/cm with the development of the disease. In the case of experimental spectra, one can notice changes in other spectral ranges. The absence of such changes in the reconstructed spectra is explained by the simplicity of the model: during the MM development, changes in the skin composition are not limited to melanin, but affect such components as triolein and ceramides [7]. The study of their influence on the Raman spectrum of skin is of interest. To determine the nature of the change in the Raman spectrum intensity with an increase in the tumor thickness, the dependences of the values in the spectrum at the Raman shift of 1575 1/cm on the tumor thickness are constructed for in vivo Raman spectra and Raman spectra obtained using the algorithm (see Fig. 5).
As can be seen from Fig. 5, in the case of spectra obtained using the algorithm, the dependence of the Raman spectrum intensity on the tumor thickness looks like exponential. With an increase in the tumor thickness, the changes become less pronounced, the intensity of the reflected light asymptotically approaches a certain limiting value -the penetration depth of the incident light into the tissue is finite and is determined by the optical properties of the tissue and the light wavelength.
For the case of in vivo Raman spectra, it is not possible to obtain the form of the dependence due to the lack of accurate data on the tumor thickness in patients whose Raman spectra were used to obtain averaged in vivo spectra. However, knowledge of the stages of the disease development allows us to conclude that the dependence is proportional to the stage.

Determination of changes in the concentrations of skin components during the development of MM
In this experiment, the geometric characteristics of the tissue and the parameters of the light source are similar to those specified in Subsection 3.1, except that the geometric structure of the optical model is presented as shown in Fig. 1, top left. The averaged in vivo Raman spectrum of normal skin is chosen as the "background" spectrum ( Fig. 2). By varying the concentrations of eight Raman active components (i.e. collagen, elastin, keratin, cell nuclei, triolein, ceramides, melanin and water), it is proposed to obtain a spectrum that is closest to the averaged in vivo Raman spectrum of skin with MM. Fig. 6 presents the results of modeling the Raman spectrum of skin with MM. The parameters at which the spectrum obtained is closest to the average in vivo Raman spectrum of skin with MM are presented in Table 2.
Comparing the reconstructed MM Raman spectrum with the average in vivo MM Raman spectrum (see Fig.  6), it can be noted that these spectra are fairly well correlated with each other, especially in the range from 1200 to 1265, from 1400 to 1440, from 1480 to 1560, from 1600 to 1700 1/cm. In the range from 1440 to 1460, from 1700 to 1740 1/cm, the obtained Raman spectrum is more similar to the Raman spectrum of normal skin. In the range from 1280 to 1300, from 1740 to 1760 1/cm, the obtained spectrum does not coincide with any of the in vivo spectra. This means that during the development of MM, changes in the composition of the skin are not limited to the components studied.
In the transition from the averaged in vivo Raman spectrum of normal skin to the averaged in vivo Raman spectrum of the MM, the concentrations of elastin, keratin, and ceramides decrease, and the concentration of collagen, triolein, nucleus, melanin, and water increases (see Table 2), that are little correlated with results in Ref. [8]. These results may be due to the fact that the composition of MM varies depending on individual characteristics of patients. The main indicator of the development of MM is an increase in the concentration of melanin, as reported in Ref. [8].

Determination of changes in the concentrations of skin components during the development of BCC
In this experiment, the geometric characteristics of the tissue and the parameters of the light source are similar to those specified in the previous subsection. The averaged in vivo Raman spectrum of normal skin is chosen as the "background" spectrum. By varying the concentrations of eight Raman active components (i.e. collagen, elastin, keratin, cell nuclei, triolein, ceramides, melanin and water), it is proposed to obtain a spectrum closest to the averaged in vivo Raman spectrum of skin with BCC. Fig. 7 presents the results of modeling the Raman spectrum of the skin with BCC. The parameters at which the spectrum obtained is closest to the averaged in vivo Raman spectrum of skin with BCC are presented in Table 2.  Comparing the reconstructed BCC Raman spectrum with the averaged in vivo BCC Raman spectrum (see Fig. 7), it can be noted that these spectra are fairly well correlated with each other, especially in the range from 1316 to 1723 1/cm. In the range from 1723 to 1799 1/cm, the resulting Raman spectrum is more similar to the normal skin Raman spectrum, and not to the BCC Raman spectrum. Paying attention to the Raman spectra of skin components (see Fig. 7), it can be noted that in the specified range the spectra of the components have no characteristic peaks. This means that during the development of BCC, changes in the composition of the skin are not limited to the components studied, the change in the in vivo Raman spectra in the specified range is due to other components.
Analyzing changes in the concentrations of the skin components, one can identify the following. In the transition from the averaged Raman spectrum of normal skin to the averaged Raman spectrum of the BCC, the concentrations of collagen, triolein, ceramides, and melanin decrease, and the concentration of water increases, as also reported in Ref. [8]. In the case of elastin, keratin and nucleus, there is a difference from results represented in Ref. [8]. Thus, in accordance with Feng X. et al. [8], the concentrations of elastin and nucleus increase by 11%, and the concentration of keratin does not change (see Table 2). However, as a result of modelling, the concentration of elastin decreases, the concentration of nucleus does not change, and the concentration of keratin increases. These results may be due to the fact that the composition of BCC varies depending on individual characteristics of patients.

Conclusions
This paper presents new approach and an algorithm for additive simulation of malignant neoplasms of human skin. Using the algorithm, Raman spectra of human skin with MM of various thicknesses (i.e. different stages) and with BCC were obtained. Obtained results were found to be consistent in general with previous studies.
Using the proposed algorithm, we can make a workaround for a problem that we know only a small part of components contained in the human skin. The process of modeling a tumor turns into a process of tweaking the component composition of the skin, simulating evolution of the tumor from healthy tissue.
The algorithm can be used not only to simulate the Raman spectra but also to study the real Raman spectra of patients and determine certain parameters of these spectra.
However, it should be noted that the Raman active component spectra used in our research were obtained using a confocal Raman microscopy [8], and the experimental skin spectra were recorded using a fiberbased Raman spectroscopy. Therefore, the results should consider this fact in a case of strong impact of the fiber in Raman spectra.

Disclosures
All authors declare that there is no conflict of interests in this paper.