2.1.

Genetic Algorithms

A standard GA consists of problem encoding, fitness function evaluation, crossover operation, and mutation operation.^{16 17} A flow chart is shown in Fig. 1. For the application of GAs to the topic under study, the iterative algorithm is performed as follows.

Figure 1

Flowchart of a GA and examples of crossover and mutation operations.

By repeating steps 2. and 3. in an iterative fashion, the GA converges toward the optimal chromosome.

The convergence time is related to both the number of parameter sets in one population (n) and number of iteration steps (N). Greater n and N will improve the likelihood of convergence but will require longer computation times. After trying different combinations of n and N for the current problem, we decided to use a small population (n=8) of randomized initial chromosomes; simultaneously, to avoid being trapped by “bad parent genes,” we allowed the algorithm to restart ten times after every 300 iterations. The best result was selected from the fitted spectra produced over these ten runs.

The advantage of GAs is that randomizing the initial starting chromosomes permits exploration of the entire solution space, and thus serves as a powerful tool to find a global, rather than a local, optimum.

2.2.

Analytical Model

Ideally, GAs may be used to optimize any computational method. However, numerical methods such as Monte Carlo simulations¹⁸ when used conventionally are very time consuming and computationally intense to be incorporated into iterative optimization loops. Recently, Spanier and coworkers have developed adaptive Monte Carlo (MC) methods^{19 20 21} that dramatically increased the rate of convergence of simulation for transport problems. However, even these improved MC algorithms will require at least several minutes to compute one spectrum over the entire visible wavelength range. In the present work, we adapted an analytical model based on an optical diffusion approximation detailed in the thesis of Keijzer.²²

2.3.

Selection of Parameters

It is necessary to decide how many skin layers should be modeled and which parameters should be included in the optimization algorithm. Since some skin parameters are strongly interdependent, it is very difficult to describe mathematically the interdependency between parameters; on the other hand, it is also impractical to include all the parameters into the fitting procedure, which will cause nonuniqueness of the solution space.

We first established a five-layered model according to Douven and Lucassen.⁵ Each layer is characterized by parameters describing the optical properties and morphology. Then we set out to reduce the number of parameters to be optimized by determining the relative impact of each parameter on the visual reflectance spectrum. For example, to determine the impact of the papillary upper dermal thickness (z2), we set all other parameters constant at their average values, and varied z2 by ±30%.³ For all z2 values, we then plotted the reflectance spectra and compared them quantitatively. Since there were only slight offsets (∼±0.7%) appearing in the 600- to 800-nm wavelength range, we decided to set z2 as a constant throughout the optimization process. By performing this procedure with each parameter, and assuming ±1% as the threshold maximum offset, we concluded the following.

1. For VRS, there is virtually no advantage in using a five-layered skin model that includes the deep blood plexus and subcutaneous fat as the fourth and fifth layers, respectively. Within reasonable parameter ranges collected from the literature,^{2 3 4 7 23 24 25 26 27 28} the parameters for the fourth and fifth layers have a negligible effect on VRS. This conclusion agrees with that previously reported in the literature. Lakmaker, Pickering, and van Gemert²⁹ reported that the maximal depth at which ectatic blood vessels contribute to the abnormal PWS skin color is 0.9 mm for a normally pigmented epidermis and 0.8 mm for a darker pigmented epidermis. Also, Niechajev and Clodius³⁰ reported histological analysis on 50 PWS, which showed that PWS blood vessels in the deeper fourth and fifth layers were virtually nonexistent. Therefore, in the present study, we use a three-layered skin model as shown in Fig. 2. (The Maple® code for solving the optical diffusion equation and the Matlab® code for calculating the wavelength-dependent visual reflectance spectrum based on given skin properties are available for download from the website in Ref. 31).

Figure 2

Three-layered skin model used in the current study, symbols corresponding to Table 1. (Schematic courtesy of Wim Verkruysse, PhD thesis, 1998.)

2. The epidermal thickness (z1) and blood volume fraction in the first layer (b1) are strongly interdependent; use of both parameters as optimization targets will cause nonuniqueness to the solution. Therefore, we chose not to vary these two parameters simultaneously. In this study, we allowed z1 to be varied and set the epidermal blood volume fraction as a constant, i.e., b1=0.2%, ⁴ to account for the few vessels protruding from the upper dermal papillary loops into the first planar layer with thickness z1.

3. Douven and Lucassen⁵ selected the average blood vessel diameter in each layer as constant average values of normal skin (8, 12, and 30 μm for capillaries, superficial plexus, and dermis, respectively). In PWS skin, the average vessel diameter in the upper dermis was reported in the range of 10 to 50 μm after 3-D reconstruction of PWS anatomy using a series of histological sections.³² A PWS blood vessel diameter of 10 to 200 μm was also reported by Svaasand et al.⁴ Since vessel diameter significantly influences light absorption³³ and laser treatment response,³⁴ we chose the average vessel diameter (specifically, the blood vessel diameter in the second layer, d2) as one of the parameters to be fitted by our algorithms. In summary, we selected nine parameters to be included in the optimization routine, and set the remaining parameters at constant values as listed in Table 1.

Table 1

2.4.

Implementation

The diffusion model was programmed as a Matlab® function. The GA was programmed as an optimization tool using Matlab®.

To focus on fitting the reflectance spectra local minima in the wavelength regions of 500 to 600 and 680 to 780 nm, higher sampling (i.e., one sample per 2.2 nm) was used in these regions, while less sampling (i.e., one sample per 4.4 nm) was used elsewhere. Similarly, when there were similar residues achieved among fitting curves resulting from ten GA runs, we selected the one with the best fit in the wavelength regions of 500 to 600 and 680 to 780 nm. The normalized residue was chosen to be the fitness function in the optimization procedure:

3.1.

Fitting Simulated Spectra

To test our method, simulated spectra obtained from forward calculation based on known parameters were used as input into the GA. Three sets of parameters were used to simulate three noticeably different spectra representative of normal and PWS human skin [Fig. 3(a)]. To simulate the same noise amplitudes as in our measurements where substantial noise was observed at relatively short and long wavelengths, white noise with amplitudes of ±0.01 and ±0.02 was added to the calculated spectrum at wavelengths of 520 to 650 nm and at all other wavelengths, respectively.

Figure 3

(a) Comparison of fitting three simulated spectra with added white noise using GA and MI methods: the fitted spectra by using GA overlap with the original input spectra very well, while the fitted spectra by using MI have much larger residues. (Retrieved skin parameters corresponding to Table 2.) (b) Residues between simulated and fitted spectra using the GA method, the ranges of residues are at the same magnitude of added white noise.

Using the GA, all nine parameters were fitted simultaneously for each spectrum. The ranges of the initial parameters were chosen as same as those listed in Table 1. The smooth solid lines in Fig. 3(a) indicate the best fits after applying the GA. Figure 3(b) shows the residues between the simulated and fitted spectra after using the GA. It can be noted that the ranges of residues are at the same magnitude as the added white noise.

For comparison, one of our authors, who has extensive experience with manual fitting of VRS measurements,³⁵ tried to fit the same three curves without any knowledge of the input parameters. He used a simple manual iteration (MI) method, which is to vary each parameter independently and then optimize them one by one. The dashed lines in Fig. 3(a) present the best fits obtained from MI. It can be noted that the residues are much larger.

Tables 2, 3, and 4 list the original inputs, fitted parameters, and their differences for the three simulated curves using the two fitting methods. Note that when using the GA, although vessel radii sometimes (as in Tables 2 and 4, line 8) have greater than 30% errors, the relative errors of all other parameters are below 20%. When using MI, even with some known parameters, the errors are much larger. Thus, we concluded that GA performed much better than MI for fitting spectra.

Table 2

Table 3

Table 4

The computation time used for a GA with 300 generations is 30 s on a Pentium® 4 computer with 2.6-GHz CPU. Therefore, it takes 5 min of computer time to fit one curve, since the algorithm is allowed to restart ten times. In comparison, MI requires 20 to 40 min of combined computer time and constant human participation.

3.2.

Fitting Measured Spectra

3.2.1.

Comparison between normal and PWS sites

Our measurements were conducted using instruments from Ocean Optics Incorporated (Dunedin, Florida): spectrometer (SD2000), integrating sphere (ISP-REF, 10-mm aperture), fiber optic variable attenuator (FVA-UV), and diffuse reflectance standards (Labsphere Incorporated, North Sutton, New Hampshire) for calibration. A notebook computer (Presario 2500, Compaq) was used for data acquisition.

Measurements were acquired from two patients with PWS skin. Patient A [Fig. 4(a)] is an Asian male with a dark red PWS on the left upper cheek. Patient B [Fig. 4(c)] is a Caucasian male with a light pink PWS on the right upper lip and cheek. All photos presented in this work were acquired with a standardized cross-polarized digital imaging system.

Figure 4

Photo and GA fitting results of VRS measurements on normal and PWS sites from: (a) and (b) Asian male with a dark red PWS, and (c) and (d) Caucasian male with a light pink PWS. The GA approach can distinguish between different skin types, it can also distinguish skin parameters between normal and PWS skin. [Note: black circles in (a) and (c) indicate the sites where measurements were taken. Retrieved parameters correspond to Table 5.]

Measurements of VRS were taken on both PWS and adjacent normal skin sites. For each measurement, the contact pressure between the integrating sphere and the skin was kept as constant as possible. Three to five measurements were taken at each site to ensure good repeatability. The averages of all measurements were calculated and used for fitting. The fitted and the averages of measured curves were plotted in Figs. 4(b) and 4(d) for patients A and B, respectively. From the numerical results listed in Table 5, it can be concluded that:

Table 5

Fitting results of VRS measurements on normal and PWS sites from two patients.
		Patient A		Patient B
		Normal	PWS	Normal	PWS
1	Residue	0.0068	0.0041	0.0081	0.007
2	μ a,694 M  (mm −1)	0.67	0.74	0.25	0.25
3	μ s,577 E  (mm −1)	42.9	40.7	45.5	35.6
4	μ s,577 D  (mm −1)	43.1	36.9	49.8	32.4
5	b2 (%)	8.8	28.2	3.4	22.2
6	b3 (%)	9.3	6.8	9.8	10.4
7	z1 (mm)	0.051	0.050	0.057	0.053
8	d2 (mm)	0.024	0.080	0.028	0.106
9	SaOxy (%)	93.2	84.6	94.6	97.5
10	Hematocrit (%)	48	45	44	46

1. (1) patient A (Asian) has a higher melanin absorption coefficient (μ _a,694 ^M ) than patient B (Caucasian), as shown in Table 5, line 2

2. (2) both patients have a considerably higher blood volume fraction in the upper dermal layer (b2) on PWS as compared to normal skin (Table 5, line 5)

3. (3) both patients have an average blood vessel diameter (d2) on PWS sites that is greater than that from normal sites (Table 5, line 8), even taking into account the ±30% error margin as calculated from fitting the simulated spectra (Tables 2 and 4, line 8).

3.2.2.

Comparison of three different PWS sites

Figure 5 shows photos and fitting results of measured reflectance spectra for three PWS patients with different lesion colors on different anatomical locations. All three patients are Caucasian. Patient 1 [Fig. 5(a)] has a light pink PWS unevenly distributed over his face; patient 2 [Fig. 5(c)] has a red PWS evenly distributed on his forearm; and patient 3 [Fig. 5(e)] has a red PWS under her left eye. The fitted and the averages of the measured curves were plotted in Figs. 5(b), 5(d), and 5(f) for patients 1, 2, and 3, respectively. Table 6 summarizes the results.

Figure 5

Fitting results of VRS measurements on three PWS patients with different lesion colors on different anatomical locations. The GA approach can distinguish skin parameters for different lesions. [Note: black circles in (a), (c), and (e) indicate the sites where measurements were taken. Retrieved skin parameters correspond to Table 6.]

Table 6

Fitting results of VRS measurements on three PWS patients with different lesion colors on different anatomical locations.
		Patient 1	Patient 2	Patient 3
1	Residue	0.0082	0.0075	0.0083
2	μ a,694 M  (mm −1)	0.27	0.23	0.33
3	μ s,577 E  (mm −1)	38.6	37.4	39.6
4	μ s,577 D  (mm −1)	35.7	31.6	29.6
5	b2 (%)	26	29	38
6	b3 (%)	3	17	28
7	z1 (mm)	0.058	0.052	0.093
8	d2 (mm)	0.052	0.096	0.094
9	SaOxy (%)	95	87	99
10	Hematocrit (%)	47	43	44

1. (1) Patients 1, 2, and 3 have the lowest, intermediate, and highest blood volume fractions, respectively. This is true for the second (b2) and third layers (b3), as in Table 6, lines 5 and 6

2. (2) The average blood vessel diameter (d2) of patients 2 and 3 are both greater than those of patient 1 (Table 6, line 8).

3. (3) Comparing the melanin absorption coefficients (μ _a,694 ^M ) from these three patients (Table 6, line 2) with patient A (Table 5, line 2), it appears that Caucasians have lower μ _a,694 ^M values than Asians.

4.1.

Possible Modeling Errors

  1. The diffusion approximation is only valid in the case of light transport when the medium is dominated by scattering, i.e., the inequality μ_a/μ_s(1−g)≪1 must be satisfied³⁶ for each of the three layers in the skin model. Figure 6 illustrates the calculation of μ_a/μ_s(1−g) corresponding to the three PWS sites of patients 1, 2, and 3 [Figs. 6(b), 6(c), and 6(d)]. For comparison, the validity condition for patient 1 normal skin site was also plotted [Fig. 6(a)]. We found that the diffusion approximation was pushed to its limit for the mathematical model to be able to fit the high blood volume fraction at the PWS sites. For example, for patient 1, the maximum value of μ_a/μ_s(1−g) in the second layer reaches 0.8 [Fig. 6(b)], and for patient 3, the maximum value in the third layer is greater than unity [Fig. 6(d)]. Therefore, the diffusion approximation is no longer valid and the absolute values of parameters obtained from the curve fitting are inaccurate and the only conclusion that can be drawn is that there tends to be more blood in patient 3’s dermal layers. Determination of exactly how much more blood is beyond the limits of the present analytical model using diffusion approximation theory.

Figure 6

Diffusion approximation validity condition for (a) patient A normal site, (b) patient 1 PWS site, (c) patient 2 PWS site, and (d) patient 3 PWS site. This condition is pushed to its limit for the case of high blood volume fraction at the PWS sites.

2. Douven and Lucassen⁵ compared Keijzer’s²² diffusion approximation model with the Monte Carlo method, which is widely accepted as the most accurate simulation method. They found that in the case of diffuse illumination, the Keijzer model slightly underpredicts the Monte Carlo results over the entire wavelength range.

3. The adoption of the layered skin structure and the assumption that each layer is a homogeneous flat medium containing vessels of the same diameter are rough approximations.

4.2.

Possible Instrumentation Limitations

  1. The finite aperture of the integrating sphere results in loss of light in the wavelength range of 650 to 800 nm, leading to a reduction in the reflectance values. Considering that oxygen saturation is the only parameter under study that can cause a steep reduction of the reflectance spectrum at these wavelengths,³⁷ this may explain why an abnormally high oxygen saturation was retrieved.

2. With our measurement system, noise in the wavelength range of 800 to 1200 nm makes it impossible to utilize the information contained in those signals. Instrumentation that gives a clearer signal at a wider spectrum would be helpful.

4.3.

Parameter Selection

Parameter selection was based on published literature, which gave quite different models for calculating skin optical properties.

1. Three different values or expressions of the baseline tissue absorption coefficient (μ_a ^T) were reported. Douven and Lucassen⁵ chose μ_a ^T as a constant to be fitted, and reported values of 0.0158 and 0.0134 mm⁻¹, which are much lower values than the 0.025 mm⁻¹ for the wavelength range of 600 to 900 nm used by Duck³⁸ and Svaasand.⁴ In our study, we adopted the expression of wavelength-dependent μ_a ^T, ^{39 40} i.e., μ_a ^T=7.84×10⁷×λ^−3.255 mm ⁻¹, to agree with other wavelength-dependent coefficients used in Keijzer’s model.^{5 22}

2. For skin melanosomes, Jacques⁴⁰ used μ_a ^melanosome =1.7×10¹¹×λ^−3.48 mm ⁻¹, while the average melanin absorption coefficient of human skin given by Meglinski and Matcher³ was μ_a ^melanin =5×10⁹×λ^−3.33 mm ⁻¹, and Svaasand et al.⁴ used the average absorption coefficient at λ=694 nm, μ_a ^melanin =μ_a,694 ^melanin (694/λ). ⁴

In consideration of the previous, it is incorrect to use the numerical values obtained from the fitting algorithm as the actual quantities of the optical properties of human skin. However, the approach of combining VRS with GA is very useful to give a qualitative comparison of some skin features between different patients/sites, or on the same patient/site during different stages of treatment. VRS may be a valuable tool to monitor treatment response when used in conjunction with a diffusion theory-based model.

We believe that the time required for curve fitting can be further reduced by: 1. using hybrid optimization methods, e.g., combining GA with some local minimization methods, such as a Simplex algorithm¹⁵; and 2. providing more specific information about the patient under study, such as skin type and PWS color, before starting the fitting procedure. Thus, the range of certain fitting parameters, such as melanin absorption and blood volume fraction, can be reduced and constrained throughout the fitting algorithms.

An as yet unsolved problem is that we find that some measured VRS, specifically those of normal Caucasian skin, are very difficult to fit. Figure 7 shows a typical fitting result from a measurement done on the forearm of a 35-year-old Caucasian male. We arbitrarily enlarged the range of each fitting parameter (as in Table 7), but the differentiation between the wavelengths 600 to 650 nm is not represented by the fitted spectrum. One possible explanation is that the skin model is too simple to include all of the absorbing and scattering factors within human skin, particularly at some wavelengths.

Figure 7

Measured and GA fitted VRS on a normal skin site of a male Caucasian. The differentiation between the wavelengths 600 to 650 nm is not represented by the fitted spectrum (retrieved skin parameters correspond to Table 7.)

Table 7