Detection of oxygen saturation in general, and minimally invasive detection specifically, is an important task. Measuring oxygen saturation is significant in many procedures, including tumor detection, cancer treatment adjustment, and ischemia monitoring during medical procedures.1 Detecting the oxygen saturation level in an efficient minimal invasive method can substantially improve treatment in such cases.
There is a variety of minimally invasive methods for estimating oxygen saturation. Most of them are based on the differences between the optical properties of oxygenated and deoxygenated hemoglobin in the visual and near-IR regions, such as absorption or scattering.2, 3
One of the possible suitable methods is photothermal spectroscopy which is sensitive to changes in the absorption spectrum. This method was investigated thoroughly and used for numerous applications.4, 5, 6 The high sensitivity of the method to surface measurements7 is extremely valuable in cancer detection, since the majority of human cancers arise from epithelial cells.8
Recent progress in the development of a flexible coherent hollow waveguide bundle in the mid-IR range made by our group and others9, 10 allows the use of such a method through a commercially available endoscope and expands the detection possibilities within body cavities. The bundle enables imaging of saturation in the examined tissue using a thermal camera, and reduces the error caused by the fiber mediation by eliminating the need of a fixed distance and orientation from the fibers.
Although there are some implementations of saturation measurement methods for internal use, none of them utilizes the advantages of the waveguide bundle and photothermal spectroscopy for oxygen saturation imaging of internal tissues.
The objective of this research was to develop a minimally invasive thermal imaging method to determine the oxygenation level of an internal tissue. In this method, the tissue is illuminated by a laser and observed by a thermal camera through the coherent hollow waveguide bundle. As a result of the photothermal effect, the tissue temperature rises. Small temperature changes can be detected because of the high resolution of the system. Since this increase depends on tissue composition, illuminating it in several wavelengths provides sufficient information for the estimation of the oxygenation level.
The objective of this work is to present an algorithm developed for assessment of oxygen saturation, and the results of implementing it on a theoretical tissue model, which has been developed for this study as well. The algorithm is examined initially on exposed skin tissue and later on internal organs.
Description of the model
A theoretical model of the problem was implemented to help design future experimental setups and develop the test procedures. The model simulates the temperature rise in the tissue as a result of the laser illumination in the same way that a thermal camera would capture it. A schematic chart of the model is shown in Fig. 1 .
There are several theoretical models that describe the optical behavior of the skin, varying on the layer division of the tissue and the method of calculating the effect of illumination.11, 12 The most detailed model, a seven-layer model of the skin, was selected to estimate the optical properties of the skin depending on depth and composition.13
The integrated model was implemented using Matlab (The Mathworks Inc., Natick, Massachusetts). All the subprograms use a 2-D model with axial symmetry assumption. Most of the parameters of the model, such as layer thickness or anisotropy, are considered as constants and were not changed between calculations. However, they can be easily changed if necessary. Other parameters, such as saturation or hemoglobin and melanin concentrations, were entered as input variables to the program.
The optical properties of each layer were calculated, taking into account material concentrations and absorption spectrum, including the baseline absorption of the tissue (the absorption of average tissue, without the main absorbers)14, 15.
The program uses a Monte-Carlo multi-layered media simulation (MCML16) to evaluate the absorption (photon density) in the skin for a given wavelength, depending on the calculated optical properties. The convolution program CONV17 was then used to calculate the absorbed energy in the tissue as a function of radial distance and depth, considering the radius of the beam and its shape.
The temperature rise of the tissue after illumination was calculated using COMSOL, a finite-element differential equation solver program (COMSOL AB, Stockholm, Sweden). COMSOL solved the heat transfer problem created using a predefined bioheat module, which is specifically designed for biological tissues. Although COMSOL has a user interface, it was integrated into the simulation using the COMSOL script to fully automate the simulation. The final temperature distribution on the surface of the skin was considered the predicted image that would have been seen by a thermal camera.
By illuminating a large region of the tissue, using a wide beam or a scanning apparatus, the saturation image of the tissue can be obtained. Otherwise, to estimate the saturation at a single point, only the temperature measurements at the point of interest are necessary. The latter option was selected for the calculations throughout the work.
The MCML program receives as input the optical properties of the layers, their thickness, and the number of photons to simulate.
Skin layer parameters used as input to the model.13
|Layer||Thickness [μm]||Refractionindex||Anisotropy||H2O concentration [%]||Bloodconcentration [%]|
|Upper bloodnet dermis||80||1.39||0.95||0.6||0.3|
|Deep bloodnet dermis||80||1.38||0.95||0.7||0.1|
The concentrations of the melanin and hemoglobin were modified between simulations, and the absorption coefficient of each layer was calculated accordingly.
The number of photons used as input for the simulation was usually 500,000 because it produced reasonable repeatability in the final results (fluctuations of less than —the resolution of the thermal camera). Since the accuracy of the model will be unknown until the experimental stage, the repeatability of results was used as a measure of accuracy.
The radial resolution and the -axis resolution were selected as 0.01 and , respectively, for the MCML program and as 0.01 and for the COMSOL program. The resolution was selected according to computational limitations of both programs and computing time limitations. The effect on the temperature calculation results was compared to temperature calculations with better resolutions, and was found to be negligible (less than the resolution of the thermal camera).
An additional difference between the models that were calculated by both programs was the tissue’s depth. Since the MCML calculations showed negligible absorption in deep layers of the tissue and their effect on external layer temperature is small, the depth of the COMSOL tissue model was reduced from to reduce calculation time and memory use.
The radius of the Gaussian laser beam used for the calculation of CONV was , and accordingly the radial size of the tissue in COMSOL was selected to be (instead of infinite as in MCML). The tissue illumination time was . The laser energy was changed according to the melanin content of the tissue to achieve tissue temperatures up to . Since the tissue is heated locally and for a very short period of time, the damage to the tissue is minimal if not negligible.18 Although the tissue is heated several times, the short and moderate heating also prevents other changes in the tissue affecting the optical and thermal parameters. The tissue’s thermal parameters were based on COMSOL’s database in correlation with existing literature. The density of the tissue and the blood was . The specific heat was . The thermal conductivity of the tissue was . Blood perfusion rate was .
Although all parameters of the model can be easily changed, only those that had a significant effect on the tissue absorption and therefore on the temperature were changed (such as melanin and hemoglobin concentration and hemoglobin oxygen saturation). These changes account for different types of tissue.
The melanin concentrations in the simulation were changed between 0 to 25% (mostly 0 to 10%). Melanin concentration of 0% is technically skin without melanin, but can also represent an internal tissue. Melanin concentration of 10% represents moderately pigmented skin.14
The minimal hemoglobin concentration in the simulation was representing a state of anemia, and the maximal value was representing a normal state for men and high concentration for women.
The hemoglobin saturation was calculated as the ratio between the concentration of oxygenated hemoglobin and the total concentration of the hemoglobin. The simulated values included the entire range of 0 to 100%.
Running the simulation for a single wavelength with a PC processor requires about . Multiple runs with changes in parameters are easy to perform and require additional per run.
Theoretical Model of the Temperature Function
The measured temperature dependence on the oxygen saturation was studied to develop the required saturation experimental measurement method.
The external temperature increase is composed of the contributions of all the heated layers in the tissue. Each layer that absorbs energy is heated, and some of the heat is transferred to the upper layers and eventually affects the external temperature. The effect of each layer depends on its depth, physical properties, and other parameters.
The total external temperature increase will be:.
The absorption of each layer is a function of its effective absorption coefficient and the available energy for absorption arriving to it. The effective absorption coefficient is the weighted sum of the absorption coefficients of the materials in each layer, according to their relative concentrations. The available energy for absorption arriving to a layer depends on the properties of the layers above it. Clearly, when the upper layers absorb most of the available energy, the temperature increase at the layer under consideration will be small, even if it has a high absorption coefficient. In a similar way, the same layer, placed below layers with small absorption coefficients, will absorb more light and will heat to a higher temperature. The effect of the upper layers depends on many variables, such as their thermal and optical properties, thickness, and other physical properties.
Therefore, the temperature function of the tissue has the form:represents the absorption of the ’th layers. The functions represent the effect of the other mentioned parameters on the absorption of the ’th layer (the absorption of the above layers and the thermal and physical properties of the ’th layer, not shown in the equation).
The function is independent of the absorption of any layer and depends only on thermal and geometrical properties. Although other optical properties, such as scattering, also depend on the wavelength, it is assumed that they are constant while dealing with a narrow wavelength range. Therefore, is independent of wavelength and can be approximated as a constant for a specific set of measurements in similar wavelengths.is independent of wavelength, except for the dependency of . Therefore, it can be approximated using a Taylor series for the variable of .
Since only is a function of wavelength, the function and its derivatives for a specific value of (in this case ) can be written as:are coefficients that need to be determined from the problem.
The function can be written in a similar way., it is possible to rewrite the temperature function shown in Eq. 2. The number of coefficients in the equation is determined by the required accuracy. For example, it can have the following form:
An example of the results of CONV is presented on Fig. 2 . The example represents normal skin tissue with 5% melanin concentration, hemoglobin, 90% saturation. The excitation was done by a cw light source for and the radius of the laser beam was . The absorption is displayed in .
Most of the energy is absorbed in the epidermis (because of the melanin) and in the upper blood net dermis (because of the hemoglobin). Deeper layers, or layers without a significant concentration of these highly absorbing materials, absorb less energy but are not insignificant. The color scale represents the absorption in units of .
The results of the COMSOL analysis for the same example can be seen in Fig. 3 .
Figures 2 and 3 show that the tissue radius that is affected by the laser illumination is lower than . As mentioned before, this value was selected as the radial size of the tissue for the COMSOL simulations. The thermal image the thermal camera will capture in this case is presented in Fig. 4 .
The time constant for temperature relaxation back to equilibrium was also calculated and is approximately .
To examine the effect of melanin in different wavelengths, a preliminary calculation was performed on large wavelength range, with low wavelength resolution . Figure 5 shows the surface temperature on the illumination point as a function of wavelength for several oxygen saturation levels (0, 25, 50, 75, and 100%) and melanin concentrations (5, 15, and 25%). Since the illumination source operates in the visual or NIR range, it is not expected to affect the temperature measurement, which is usually performed in the range.
The figure shows that the effect of the saturation on the temperature decreases as the melanin concentration increases. The differences created by saturation changes in tissues with higher concentration of melanin are smaller and difficult to detect, and therefore saturation evaluation will be limited by the melanin concentration. When the method is applied to internal organs, the melanin concentration is irrelevant.
The saturation influences the temperature at wavelengths where there is a significant difference between the oxygenated and deoxygenated hemoglobin absorption spectrum. This influence will be negligible if the melanin absorption in those wavelengths is much higher than that of the hemoglobin. Therefore, the temperature differences created by the saturation at are smaller than those created at , in spite of the larger differences in the hemoglobin absorption in that wavelength (even after increasing the excitation energy to match the temperature level to the temperature at ).
To develop a method as suitable as possible to high melanin concentration, the most suitable wavelength range is . For internal organ applications where melanin is not present, one might consider the near-IR range, where the absorption differences between oxygenated and deoxygenates hemoglobin are higher.
The focus of the rest of this work is on the wavelength range between 410 and . However, the goal of this work is to develop a generic method that would also apply to other wavelength ranges and applications.
Comparison between Calculations and the Theoretical Model
To derive the calculated temperature dependence on wavelength, the temperature function was evaluated again, this time from calculations of the simulation on various types of skin tissues.
With that aim, the dependence of the surface temperature function on various variables was calculated.
Figure 6 shows the temperature as a function of wavelength in skin tissue with 5% melanin concentration and hemoglobin concentration with 50% saturation. The power used was for , with and without the absorption of water (the absorption coefficient is set to zero).
The figure shows that in the relevant wavelength range, the water’s absorption hardly affects the temperature because of its relatively small absorption coefficient. Therefore, the temperature function is practically independent of the water concentration.
To evaluate the temperature’s dependence on hemoglobin concentration, the simulation was performed on the same tissue model with different hemoglobin concentration in physiological levels with several melanin concentrations. The results are shown in Fig. 7 and show linear dependence, meaning that the temperature function is of the form:and are coefficients independent of .
Similarly, different saturation levels were simulated to evaluate the saturation’s effect on temperature. Figure 8 shows the temperature dependence of the oxygenation for 5% melanin concentration and hemoglobin concentration.
The results show that the temperature function is of the form:is the effective absorption of the hemoglobin, depending on the saturation: and are the deoxygenated and oxygenated hemoglobin absorption coefficients. and are coefficients independent of , , and hemoglobin absorption.
Figure 5 shows that the melanin concentration is the main contributor to the temperature increase. The simulations shown in Fig. 9 show that the melanin adds a nonlinear contribution to the temperature rise. Figure 7 shows that the changes in the melanin concentration also have an effect on the hemoglobin contribution.
Assuming the concentration and absorption coefficient of the melanin affect the absorption in the same manner (since the absorption is proportional to their multiplication), small changes in the absorption coefficient will have a linear effect on the temperature, for constant hemoglobin absorption, and negligible effect on the hemoglobin’s contribution.
Since most of the tissue is considered baseline tissue, the concentration of that tissue will not change significantly due to changes in other material concentrations. The main reason for change in the baseline absorption is the dependence of the absorption on the wavelength. The baseline absorption is at and at , a ratio of 1.5. Figure 10 shows the same calculations as in Fig. 7, with doubled baseline absorption coefficient. The hemoglobin contribution was hardly changed and the total temperature contribution of the baseline tissue depends linearly on the melanin concentration.
The re-evaluated temperature function is therefore of the form:and are the melanin absorption and concentration, respectively. and are the baseline tissue absorption and concentration, respectively. and are functions of melanin and baseline material concentrations and absorption coefficients.
However, for a narrow wavelength range, the temperature function can be approximated to the following simple equation:is a constant and is a bilinear function of the melanin and baseline absorptions.
Following the previous calculations, the following assumptions can be made:
1. The calculations of the effective absorption of the different layers do not need to take into account the water content of the tissue, and may include only melanin, hemoglobin, and baseline absorption.
3. By comparing Eqs. 2, 12, it is seen that is approximately the hemoglobin absorption (the hemoglobin is the dominant absorber in the upper blood net dermis, the forth layer). Any effect of the baseline absorption in the layer can be corrected by slightly changing the bilinear function .
4. Since the concentrations are unknown and are multiplied by other unknown coefficients, the equation can be written in the following simple form:8 (because the melanin layer is indeed above the baseline layer and the stratum corneum absorbance is negligible).
The temperature function has a known form depending on a finite number of unknown coefficients ( , , , , , and ) and a set of well-known variables (the absorption coefficients ). The unknown coefficients are constant for a specific set of measurements (although they can change between people and measurement locations), and therefore by measuring, or in this case calculating, the temperature in a sufficient number of wavelengths, there will be sufficient data to extract the unknowns, specifically the saturation.
Although the temperature function is known, direct solution of the equation system is not recommended due to its high vulnerability to measurement errors and inaccuracies in the theory. Even though this method might work for the results of the simulation, it might produce poor results on real experiments, unless the function is calibrated and evaluated again. This would require measurements for a wide set of material concentrations and different skin types, and could be done at a later stage if necessary.
Therefore, a more generic method is preferred. A generic method could also be applied to internal organs with minimal adaptation.
The selected method was curve fitting between the measured temperature curve and the theoretical temperature function, with the coefficients and the saturation value as unknowns. The algorithm inputs are the estimated temperature function and a set of measurements. Its output is the specific set of unknowns that provide the best fit to the temperature equation. An example for a suitable algorithm is Matlab’s curve-fitting toolbox, which was used throughout this work.
The accuracy of curve-fitting algorithms improves as the number of measurement points increases (in this case, the number of temperature measurements). However, additional measurements require more time. Therefore, the number of measurements should be a compromise between time and accuracy constraints. In this work, each set of calculations included nine different wavelengths. As mentioned, the selected wavelength range for the saturation evaluation was , because in that range the hemoglobin has a higher absorption coefficient compared to the melanin. The calculations were performed for the following wavelengths: 410, 414, 418, 422, 426, 430, 434, 438, and .
The probability of finding a solution (finding a set of values that would cause convergence of the algorithm error to zero) depends on the number of unknowns and on the initial guess entered into the program. Therefore, reducing the number of unknowns and developing a good initial guess algorithm are essential to finding the saturation accurately.
To decrease the number of unknowns, the input function to the algorithm was selected as the difference between the temperature function at the measured wavelength and its value at , which eliminates the need to find . The new function, marked as , is therefore:(10%) to reduce the number of unknowns. Since the algorithm searches for the best fit to the actual function, it will rectify the error caused by changing the weight of the papillary dermis baseline absorption accordingly.
Calculating the temperature function’s derivative might have low accuracy, depending on the wavelength resolution of the temperature measurements and their accuracy. However, using this equation in the curve-fitting algorithm is less sensitive to the initial guess, and therefore it can be used as the initial guess algorithm. Figures 11 and 12 show the results of the algorithm when applied on tissue temperature calculations for different concentrations of melanin. The initial guess of the saturation was 50%. The estimated saturation is displayed as a function of the true saturation that was used for the simulation. The black line indicates the ideal result where the estimated saturation is equal to the true saturation.
Table 2 shows the root mean square (rms) of the error of the saturation estimation results of the final algorithm.
Error rms of saturation for different skin compositions.
Skin tissue without melanin—internal tissue
As mentioned, the described method can also be applied on internal tissues, using a commercially available endoscope. To evaluate the accuracy of the method when applied on such a tissue, it was tested on the calculated skin tissue model, without melanin. This model might not accurately describe any internal tissue, but as a generic sample for a generic method, it should be sufficient for initial evaluation.
The temperature function is obviously different in this case. The melanin, which was very dominant in the skin tissue, does not exist in the internal tissue modeled. Therefore, the absorption of the epidermal layer will now be approximately equal to the baseline absorption. In addition, the strong absorption of the melanin masked the effect of lower layers, and without it, they become more significant. Therefore, the absorption of the deep blood net dermis is stronger and affects the temperature increase. The absorption of this layer is affected by the above absorption of the hemoglobin and baseline tissue. According to the previously described approximations, that would be manifested in the temperature function seen here.
Figure 15 shows the result of the saturation initial guess algorithm.
The results of the initial guess algorithm are poor, but they were sufficient for receiving good results in the final algorithm, and were therefore chosen.
Figure 16 shows the result of the final algorithm.
Table 3 shows the root mean square of the algorithm error.
Error rms of saturation for skin tissue without melanin (internal tissue).
An experimental setup implementing the method was set up in the laboratory. The models were illuminated by a continuous-wave tunable titanium-sapphire laser (3900S, Spectra Physics, Newport, Mountain View, California) pumped by a frequency-doubled continuous-wave neodymium/yttrium-vanadate (Nd:YVO4) laser (Millennia Vs, Spectra Physics), and were imaged by a thermal imaging camera (Thermovision A40, FLIR Systems, Wilsonville, Oregon) operating in the spectral range.
The imaging resolution is the camera’s spatial resolution, which is . The thermal sensitivity (affected by the spectral resolution and other factors) is , which is of the same order as the simulation errors and therefore reasonable. The camera acquires images at , faster than the heating (which is on the order of a minute), and is therefore more than adequate for this application.
The relative intensity of the laser at every wavelength was measured and used to normalize the results to an equal intensity level. The measurements were performed in several wavelengths between 720 and .
Phantoms were prepared by adding various absorbers to agar solution (Gibco, Carlsbad, California). The absorbers were used to simulate materials with different absorption coefficients, representing the oxygenated and deoxygenated hemoglobin and the melanin. Since the method is not limited to estimating the saturation of a specific material, the absorption spectra of the materials used does not have to be similar to the spectra of the materials they represent.
The hemoglobin was represented by ICG (Cardiogreen, Fluka, Milwaukee, Wisconsin) and methylene blue (Sigma-Aldrich). The absorption spectra of ICG and methylene blue solutions were measured by a spectrophotometer before each set of measurements to get accurate input for the algorithms (as the ICG spectrum changes in different concentrations). The melanin was silmulated with ink solution (Indian ink, Talens, Netherlands). The absorption coefficients of all materials—methylene blue, ICG, and the ink—can be seen in Fig. 17 .
The theoretical simulation used does not consider model cooling by external air, and therefore in the experiments, the heating was slower and the temperature does not exceed a certain limit. This limit temperature is reached when the model is in thermal equilibrium with the surroundings. The temperature increase is considered as the difference between the limit temperature and the initial temperature.
Since the heating is slow, the limit temperature is reached after several minutes. An algorithm decreasing the required measurement time was developed to prevent inaccuracies resulting from the long measurement (changes in external temperature, calibration drift, and model bleaching). This algorithm estimates the initial temperature as the average temperature before the illumination, and the heating as exponential using curve-fitting methods. Using the algorithm shortens the required measuring time and reduces the effect of noise on the temperature estimation. Figure 18 shows the original temperature measurement (solid line) and the estimation (dashed line) for the 100% ICG model during illumination. In this case the initial and limit temperatures were estimated at 293.57 and , respectively.
The temperature increase for each wavelength used for the illumination was calculated in this manner and used to estimate the model’s saturation. The saturation is defined as the ratio of methylene blue to (methylene blue + ICG).
Initially, single layer phantoms were prepared and measured in five wavelengths (720, 740, 760, 780, and ). The initial saturation guess was 50% and the temperature function used was:19 . The solid lines represent , 0, and deviations from accurate estimation.
The accuracy of the algorithm was reasonable (8.5% rms) and therefore more complicated phantoms were created. These phantoms were composed of two layers: an absorbing upper layer simulating the epidermis that contains ink, and the lower layer simulating the blood net layers containing the ICG and methylene blue.
The temperature function used was:) were used. To improve the estimation, the results of the saturation estimation were used for additional runs of the algorithm to reach better convergence. In this case, instead of 50% initial guess, ten different initial guesses were used and the final results were averaged to achieve to final estimation. The results are shown in Fig. 20 .
Results are in agreement with the real saturation values, except for the 25% model. The large error is due to malfunction in the temperature measurement in the wavelength. The rms of the error is 15% (9.2% without the problematic model).
The theoretical simulations show that for tissue with small to medium concentrations of melanin, the rms of the error is 5% to 9%. This is less accurate than available commercial methods, but it is still remarkably good for a preliminary method. In addition, the results show good feasibility for developing this method to the required accuracies. This is further emphasized by the fact that these results are achieved without any calibration and are all over the saturation range.
Examining the results for different concentrations of melanin and hemoglobin shows that the accuracy for 5 and 7.5% melanin is slightly better than for 2.5% or 10% concentrations. The decrease in accuracy for very low melanin concentration can be explained by the smaller absorption of the melanin compared to the baseline absorption, resulting in an intermediate state between the different models. This problem can be solved by using a more complex model or by choosing a model according to skin color. This option might not be suitable if the method is meant to be used without any calibration.
High melanin concentrations cause substantial absorption by the epidermal layer, and therefore a decrease in the energy available to absorb in deeper layers and smaller absorption by the hemoglobin. The increase in the ratio between the melanin absorption and the hemoglobin absorption therefore causes masking of the hemoglobin absorption by that of the melanin, and decreases the accuracy of the algorithm.
Due to the same reason, very low hemoglobin concentration can also be masked by the melanin, and the accuracy is expected to decrease.
Obviously, some of the perturbations of the accuracy, seen in the table, are caused by other inaccuracies such as the inaccuracies of the tissue model, errors in the calculations, and errors of the method itself.
Experiments on phantom models show that the method can be easily implemented and that the algorithm can be adapted to estimate relative concentrations of other materials. The algorithm gives accurate results as long as the temperature is measured properly and the developed noise reduction algorithm is used.
The number of wavelengths used is a function of the complexity of the model, from five in the simple phantom models to nine in the complex theoretical simulation models.
Further improvements of the method and the model to increase accuracy are not implemented at this stage. In the future, the method will be tested in in-vitro and in-vivo experiments and improved accordingly.
Future research will also consider the effect of the errors added by the fiber imaging bundles, the thermal camera, and other elements of the system, and possible methods for decreasing them.
The nature of the method makes it a good candidate for saturation imaging in internal cavities. For example, it could be used to detect early stages of colon cancer during a simple colonoscopy. External applications are also possible, such as detection of cancerous regions in the breast, or any other tumor close to the surface.
Calibrating the method (for example, in in-vivo experiments) will allow measurement of the concentrations of the materials and not only relative concentrations, such as saturation estimation.
A new minimally invasive method for oxygen saturation measurement is presented. This method can be used to evaluate the oxygen saturation of a tissue, external or internal. The method does not require calibration or knowledge of the specific structure of the tissue, only of the main absorbers in it. The method enables imaging of the examined tissue, thus allowing detection of ischemic tissues or cancerous regions, and can improve early detection rates and help adjust cancer treatments. The accuracy of the method is good, considering the early stage of research and the mentioned advantages.
The authors would like to thank the US-Israel Bi-National Foundation for its generous support through Grant No. 2007382.