Multilayer diffusion modeling and Coherent anti‐Stokes Raman scattering microscopy for spatially resolved water diffusion measurements in human skin

In this work using Coherent anti‐Stokes Raman Scattering microscopy, it was possible to directly measure the time dependent, spatially resolved change in concentration of water (D2O) in intact skin tissue with a spatial resolution of under 1 μm, and combined with a multilayer diffusion model, diffusion coefficients at different depths in the tissue were extracted. The results show that the diffusion varies at different layers throughout the Stratum Corneum (SC), indicating that the SC is not a homogeneous barrier but a complicated heterogeneous structure. Interestingly, averaging over the diffusion at the different depths and samples gave a relatively constant value of 0.047 ± 0.01 μm2/second. Treating the skin with acetone or tape stripping led to an increased diffusion coefficient of 0.064 ± 0.02 μm2/second and 0.079 ± 0.03 μm2/second, respectively. The combined method and model presented here shows potential for wide applications for measuring spatially resolved diffusion of different substances in a variety of different samples.


| INTRODUCTION
Recently, there has been a growing interest and realization of the need to understand cells, cell plasticity and function in the context of their surrounding intact tissue. This has led to advancements in fields such as spatially resolved transcriptions and proteomics, which deliver unique information on the complex spatial organization of tissues [1][2][3][4][5]. These results motivate the need to also study biophysical properties of samples with high spatial resolution to be able to complement the information coming from, for example, spatial transcriptions. An example of this is the spatially resolved diffusion of water or other substances in cells and tissue. Diffusion and transport of water and other materials are vital for all biological processes and take place within tissues with complex and varying spatial compositions. For example, human skin, which must maintain the homeostasis of the human body as well as providing barrier properties toward outer influences such as chemicals and foreign pathogens, consists of several different layers the outermost is being the epidermis. The epidermis also consist of several layers and the outermost layer is the stratum corneum (SC) (usually from 10 to 20 μm thick), which is said to hold the main barrier function of the skin [6][7][8]. The keratin filled corneocytes (flattened non-nucleated cells) that make up the SC are described by the well-established brick and mortar model [9]. The corneocytes (ie "bricks") are surrounded by a protein envelope comprised of crosslinked proteins covalently bound to a layer of ceramides (mortar) [8]. Recently, there has been a growing interest and realization of the need to understand cells.
Therefore understanding the barrier function as well as the transport through human skin with spatial resolution is vital in relation to environmental protection of the skin [10], skin diseases [11] and transdermal drug delivery [12][13][14]. A very fundamental process is the hydration of the SC and is said to influence several functions such as the mechanical toughness, barrier functions and regulation of enzyme activity [15][16][17][18]. The study of the hydration of the SC is therefore of great importance. Previously, the permeation of water through the SC has been studied by measuring the transepidermal water loss (TEWL) [19], sorption experiments [20] and nuclear magnetic resonance [21]. This methods however, do not provide information regarding the depth-dependent permeation. In order to investigate the depth-dependency TEWL has been performed where the SC was removed by tape-stripping [22][23][24], but this approach is destructive to the sample. In recent years, the hydration of the SC in human skin has been studied using Raman spectroscopy [17,18,25], and confocal Raman microspectrometer [15], which are noninvasive. However, these techniques offer a poor spatial and temporal resolution.
In order to investigate the transport of water through human skin we propose to use the technique of Coherent anti-Stokes Raman Scattering (CARS) microscopy [26][27][28], which allows the water in the sample to be imaged at near-physiological conditions with a resolution down from 200 to 300 nm and where optical sectioning of the sample is also possible [29,30]. CARS microscopy enables direct imaging of the penetration of water in the skin by visualizing the vibrational modes of specific chemical bonds of the molecules of interest, without the use of labeling [31]. CARS microscopy has previously been used in skin related studies to the structure of the epidermis [28,32,33] and penetration and distribution of different compounds in the skin [34,35]. In the present work, to visualize and measure the penetration of H 2 O through human skin, D 2 O was added to the surface of excised skin samples, and the penetration was imaged over time using CARS microscopy. D 2 O was used as an alternative to H 2 O as this enabled visualization of the penetrating water distinctly from the water already present in the sample. Due to the similar structure and size of D 2 O compared H 2 O, D 2 O is an excellent analog for diffusion measurements of H 2 O in human skin. Data was collected as a time series of X-Z scans, with the x axis parallel to the SC surface and the z axis perpendicular. Figure 1A shows a schematic of the sample setup and Figure 1B shows a sketch of the D 2 O penetration into the skin over time.
In order to understand transdermal delivery multiple mathematical models of diffusion in skin have been proposed with varying complexities [36], from simple macroscopic 1 dimensional homogeneous layered models [37] to more complex microscopic 2D and 3D brick and mortar models with varying complex 3D structural representations of the skin [38,39].
In this work, we present diffusion models and coefficients of transport of D 2 O through the SC obtained by using CARS microscopy. Three different mathematical models are tested in order to describe the diffusion of water through the SC. One where it is assumed that SC consists of a single layer with a single diffusion coefficient and another where SC is assumed to consist of three separate layers with three different diffusion coefficients and finally a model with multiple layers and diffusion coefficients ( Figure 1C). The goal of this being to determine if the SC can be considered as a single homogeneous membrane or if it is subdivided in multiple layers with different diffusion/barrier properties. The measurements were carried out on untreated excised human skin, skin treated with acetone and tape-stripping.

| MATHEMATICAL MODEL FORMULATION
In the following, the three models used in this work are presented, see Figure 1C. The mathematical models of the SC assume homogeneity in the x-y plane along the surface of the SC. The model is thus effectively with one spatial dimension, z, perpendicular to the surface. We place the origin, z = 0, at the outer surface of the SC, and let z grow as one moves inside. The outer surface of the SC is defined as where the diffusion in the SC becomes different than that in the donor region. Finding this height is important to be able to fit the data correctly. We model the concentration C = C(z, t) of D 2 O as function of z and time t in this infinitely broad SC region of thickness h. The region can be subdivided into several layers, inside which the dynamics of C follows the diffusion equation with diffusion constant D i , which can vary between the layers labeled by i, but is constant inside a layer. In all our models, we have Dirichlet boundary conditions at z = 0 and z = h. At z = 0, we assume that the concentration C follow a given concentration C b where A, k 1 and k 2 are found by a least squares fit to the concentration profile at z = 0. At z = h, we assume that D 2 O can quickly diffuse away, leading to C(h, t) = 0. What happens in the SC region follows three different models: the one layer model, the three layer model and the multilayer model. The one layer model: Here, there is only one homogeneous layer with the same value for the diffusion constant D 1 everywhere. The three layer model: Here, the SC region is divided into three layers of different diffusion constants, with boundaries at 0 < z 1 < z 2 < h. The extra internal boundary conditions are continuity of concentration: Þ for i = 1, 2 and conservation of D 2 O: The multilayer model: Here, the diffusion constant D i of each 1 μm layer is fitted individually. This is achieved by fitting the concentration profile at its boundaries z i and z i+1 to sums of two exponentials, that is, C(z i , t) is least squares fitted by where A i , k i1 and k i2 are the fitting parameters. The diffusion constant D i is then found by fitting the solution of three different models used to fit the data. All of the models assume an infinite dose experiment with a perfect sink condition. They also all have an outer SC layer, with a height z OSC , in which the diffusion is considered to be infinitely fast as in the donor compartment. The single-layer model describes the SC as a single layer with the height: h and one diffusion coefficient, D 1 . The three layer model describes the SC as a three layer system with the layer heights: z 1 -z 3 and three diffusion coefficients, D 1 to D 3 . Lastly, the multilayer model describes the SC as multilayer system with n layers each with height z n and diffusion coefficients D n the diffusion problem with Dirichlet boundary conditions at the layer boundary, that is, C(z i , t) = C i (t), such that the concentration at the midpoint of the layer z = (z i + z i+1 )/2 minimizes the squared deviation from the measured values there. Note that, this model does not ensure conservation of D 2 O between the layers, but instead aims at finding local diffusion constants with relatively small computational effort. Thus, it is not a physically realistic model in the same sense as the one and three layer models are.
The function C(z, t) at an arbitrary depth z is found numerically by transforming to Laplace space via where the variable name u is used to indicate that the function has been Laplace transformed. In Laplace space, the solution to the diffusion equation in region i is where the integration constants E i and F i are found via the boundary conditions by solving the corresponding set of linear equations. The concentration profile at time t is obtained by inverting the Laplace transform numerically using the Talbot algorithm [40]. As the CARS signal is proportional to the number of emitters/ concentration squared, the data is fitted to the concentration profile squared. While this is not a perfect approximation due to the mixed term between the resonant and nonresonant term it is considered to be a reasonable approximation due to the large resonant signal from the D 2 O. Furthermore, a deviation from the squared relation between the CARS signal and emitters would result in a slight increase in the measured diffusion coefficients [41].
3 | RESULTS Figure 2 shows CARS microscopy images of cryo-frozen and sectioned human epidermis where the vibrational mode at the CH 2 -stretch (lipids, 2850 cm À1 ) and at the OH-stretch of water (3356 cm À1 ) are visualized. The intercellular lipids of the SC can clearly be seen in the left panel, which visualizes the CH 2 -stretch, as the layered structure at the top of the image. Under the SC the rest of the epidermis can be seen. The dermis can be seen as the more structured part at the lower part of the image. The middle panel, depicting the OH-stretch, also shows a layered structure in the SC. However, closer examination (last panel) shows that these layers are offset relative to the intercellular lipids; clearly showing the layered brick and mortar structure of the SC where the layers of intercellular lipids (red) surround the more hydrophilic corneocytes (blue).
To measure the penetration of H 2 O through human skin, D 2 O was added to the surface of excised skin from three donors, with three different treatments of the skin, and the penetration was imaged over time using CARS microscopy. The D 2 O was added using a syringe pump attached to a needle placed under the skin. Imaging was performed close to the needle tip but without having the needle in the field of view (see Figure 1A). For the measurements the CARS lasers were tuned to the OD-stretch at 2380 cm À1 . Data was collected as a time series of X-Z scans, with the x axis parallel to the SC surface and the z axis perpendicular. The images were 512 Â 340 pixels and the scan speed was 1.3 seconds per image. The images of the D 2 O penetration were recorded for the different treatments using samples from the three donors resulting in 9 to 15 independent measurements per treatment per donor. of D 2 O into the skin before application, upon application and after 5 minutes. Before addition of D 2 O a none resonant background signal from the tissue is present. An increasing intensity (orange signal) clearly shows the addition of D 2 O, and that D 2 O is penetrating into the SC over time. Figure S1 show the same for skin samples treated with acetone and tape-stripped samples.
Intensity data from the images was collected at different depths (every 1 μm) as a function of time from the addition of D 2 O and the nonresonance background was subtracted. Examples of this can be seen in Figure 4A,B. The intensity data reflects the concentration of D 2 O at different depths as a function of time.
To further interpret the data and investigate the diffusion coefficients at different depths in the SC, the penetration data was then fitted to the three different macroscopic diffusion models, the single-layer, the threelayer and the multilayer model. The goal of this being to determine if the SC can be considered as a single homogeneous membrane or if it is subdivided in multiple layers with different diffusion/barrier properties. A diagram of the three models can be seen in Figure 1C.
For all models and samples a layer of 1 À 3 μm at the surface of the SC was found with a very fast diffusion. This fast diffusion was, in the context of the temporal resolution of the present measurements instantaneous. This The penetration of D 2 O into the SC at different depths (black) and fits made using the single-layer model (magenta) and the three layer model (green). It is clear that these models cannot successfully fit the data. (B) The penetration of D 2 O into the SC at different depths (black) and fits made using the multilayer model (red). This model is seen to successfully fit the data. In this case, all the fits have a RMSD of <0.08. The data is representative of more than 120 individual measurement areas is believed to be due to the fact that the outermost layer of the SC is preparing for desquamation and is therefore not tightly attached and free diffusion of D 2 O between these layers is possible. Similar behavior has been seen for diffusion of lipid nanoparticles in the outer layer of the SC [42].
The single layer model and three layer model could not satisfactorily fit the diffusion data at the different depths through the SC (see Figure 4A). The single layer model was typically used to fit data from 3 to 9 different depths spaced throughout the samples. Adding more data from more than three layers did not improve the fitting. The three layer model in it nature was used to fit data from three depths spaced throughout the sample. The model could simultaneously fit data from 1 to 3 depths in each of the three layers. Adding more data to the fitting process did not improve the results. It can be seen from Figure 3 that the single layer model could not fit all three layers well. In fact, both the layer at 0.5 μm and 3.5 μm are not fitted well. This is seen by the inability of the fit to follow the rise in intensity. Similar behavior is seen for the three layer model at 0.5 μm. In general, these models could not fit the data at multiple depths simultaneously.
The multilayer model, on the other hand, could fit the data well. An example of this is seen in Figure 4B where the diffusion data is shown to be fitted at multiple depths with a normalized root-mean-square deviation (NRMSD) of <0.04, which was typical three to eight times better than for the other models. The diffusion coefficients for the different layers can be seen in the figure legend. Where the values range from 0.09 to 0.04 μm 2 /second in the sample, showing clear regions with varying diffusion at different depths. This suggests that the diffusion and/or barrier properties of the SC are not constant throughout the SC. This was typical for all samples measured.
All the data from the untreated, tape-stripped and acetone treated skin samples where fitted using the multilayer model. The results of the averaged diffusion coefficients for the different treatments as a function of depth in the SC are shown in Figure 5.
Looking at the results for the normal skin (black bars), it can be seen that the outermost 2 μm have a fast diffusion which cannot be measured with the temporal resolution used to record the data. Interestingly, the averaged diffusion coefficients at different depths below 2 μm in the SC seem to be relatively constant. After averaging over the different depths below 2 μm a somewhat constant diffusion coefficient of 0.047 ± 0.01 μm 2 /second is found (Table 1). This value is comparable to the 0.216 ± 0.114 μm 2 /second found using TEWL measurements [23] and the 0.03 to 0.08 μm 2 /second found by Blank et al [43].
So while the fitting of the individual measurements clearly showed a variation in the diffusion and barrier properties at different depths, indicating that the SC is not a homogeneous barrier, on average this averages out to the relatively constant value of 0.047 ± 0.01 μm 2 /second. The difference in diffusion between the different donors can be seen in Figure 6, where the averaged diffusion coefficient through the SC along with the SD of the diffusion coefficients measured at different depths is presented. Although donor 3 is seen to have a somewhat lower average diffusion coefficient of 0.036 ± 0.013 μm 2 /s no statistically significant difference between the donors was found.
To investigate the effect of solvent induced cutaneous barrier disruption, acetone treated SC samples were investigated (magenta bars Figure 5). The outermost layer of the SC, with the fast diffusion, was found to be on average 1 μm deeper than for the nontreated skin, while the overall behavior of the measured diffusion coefficients at different depths was very similar to the nontreated skin. However the SD was three to four times larger for the individual depths. The averaged diffusion coefficient for the acetone treated skin was 0.064 ± 0.024 μm 2 /s, which was found to be statistically significantly different from the untreated skin (P « .001). The results suggest that only the outermost layer of the SC was strongly affected by the acetone and that no substantial changes in the SC barrier function were induced. This conclusion is supported by the work of Barba et al [20], who showed that acetone treatment of porcine SC did not lead to substantial extraction of SC lipids and in fact F I G U R E 5 Bar chart showing the diffusion coefficient at different depths in normal skin (black), skin treated with acetone (magenta) and skin treated by tape-stripping (green) according to the multilayer model. The dashed lines at 0.5, 1.5 and 2.5 μm indicate that the diffusion was faster than the current technique could measure. The data is representative of more than 42 individual measurement areas for each different sample types resulted in a similar or slightly improved barrier function with regard to water transport.
In order to mechanically induce a disruption of the SC barrier function measurements were performed on tape-striped SC samples (green bars Figure 5). In general the diffusion in the tape-stripped samples was slightly higher than the nontreated samples (averaged diffusion coefficient 0.079 ± 0.031 μm 2 /s), however with considerably higher SD. The average diffusion coefficient was found to be statistically significantly different from the untreated skin (P « .001) This could be due to the tapestriping process damaging the lower layers of the SC creating areas/pores with low cohesion and therefore locally damaged barrier properties. The tape-stripping procedure is expected to remove between 2 and 3 μm of the SC [22]. This has been shown to increase the TEWL in human skin by a factor of 1.4. As the TEWL is directly proportional to the diffusion in the SC, it is interesting to compare the ratio of the diffusion coefficient of the tapestripped and normal skin (0.079/0.047μm 2 /s = 1.7) to the 1.4 from the TEWL measurements. Considering that the measurements are from different donors and made using different techniques, the overall size and order of magnitude of the effect of the tape-stripping is very much the same for both measurement techniques.
Interestingly, tape-stripping combined with TEWL measurements have shown that the barrier function of the SC with regard to TEWL was constant and independent of thickness and that removing a certain fraction of the thickness of the SC had the same effect on the TEWL, independent of SC thickness [44]. This could mean that the barrier is not at the same absolute depth within the SC, thus explaining why the local differences in diffusion coefficients average out for the CARS measurements at different depths.

| DISCUSSION
In this work, a method and model to directly measure the diffusion of water (D 2 O) or other molecules in human skin tissue at different depths using CARS microscopy was successfully implemented. The intensity of the CARS signal depends on the absorption of the skin and therefore the depth into the skin. Consequently, the method can give quantitative information on the temporal behavior but not the absolute concentration of D 2 O at different depths.
CARS microscopy combined with the multilayer model enabled spatial resolution of the transported molecules and extraction of diffusion coefficients down to a resolution of under 1 μm in the sample. Thus advancing the understanding of biophysical properties of the tissue with high spatial resolution in intact tissue. A major advantage of this method compared to other methods, used to measure diffusion in cells and tissue, such as the fluorescence based Fluorescence correlation spectroscopy [45] and Raster image correlation spectroscopy [46] is that CARS is a label free technique. Therefore, it does not need an added fluorophore and can directly measure the diffusion of the molecules in question.
The spatial resolution of the technique is limited by the resolution of the microscope in the axial direction [47] and also by the depth into the sample as the point spread function will deteriorate deeper in the sample. The temporal resolution is decided by the sampling time (the time it takes to take an image) and determines which range of diffusion coefficients can be measured. By increasing the sampling rate faster diffusion processes can be measured. In practice, the sampling rate would be the limiting factor for the diffusion coefficients which T A B L E 1 Showing the average diffusion coefficients (μm 2 /second) over different depths, SD and P-value in normal skin, skin treated with acetone and tape-stripped skin

Skin treatment
Average diffusion coefficient (μm 2 /s) SD (μm 2 /s) P-value could be measured. For example, a sampling rate of 1 to 10 images per second would allow measuring diffusion coefficients of up to 1 μm 2 /s, which is comparable to the diffusion of protein aggregates and transmembrane proteins. To increase the sampling rate line or circle scans could be used instead of imaging. This would decrease the spatial information but would allow sampling rates of 1 to 8 kHz with a resonant scanner. This could measure diffusion coefficients on the order of 10s of μm 2 /s, corresponding to typical values of cytoplasmic proteins [46]. The label free visualization of CARS microscopy would support measurement of a wide range of substances such as drugs and metabolites [34,35,47]. Interestingly, it was also found, that it was only possible to fit the data satisfactorily to the multilayer model, and not the one-and three-layer models. This failure of the one-and threelayer model to fit the data underlines the need for a model, which can take the spatial variation of the diffusion properties into account. The multilayer model presented here, while not being a fully physical model, could quickly and satisfactorily fit the spatial resolved data and extract diffusion coefficients at different depths. The spatial resolution of the method revealed that in the outermost 2 to 3 μm of the SC diffusion was very fast approaching free diffusion. Indicating that the outermost layer of the SC is loosely packed and water can freely diffuse in this area. As mentioned, it was only possible to fit the data satisfactorily to the multilayer model. This demonstrates that the SC is not a homogeneous material and in fact the diffusion was found to vary at different layers throughout the SC. This is in fact quite reasonable considering that the SC is a layered structure consisting of layers of keratin rich corneocytes surrounded by the lipid matrix [8]. The results imply that the SC is clearly not a homogeneous barrier but a complicated heterogeneous structure and should be treated as such when designing and modeling transdermal drug delivery.
It is also worth to consider that the keratin in the SC can swell upon hydration and that the diffusion rate in the SC could be changed by this. As the measurements presented here, take on the order of 6 minutes a swelling of <10% is expected [48]. With regards to the diffusion rate, this is predicted to be relatively constant within the time frame of the experiment, as models have shown that the diffusion is quite constant for tissue hydration up to 0.6 g/cm 2 and it is not expected that the hydration will reach above this within the time frame of the experiment [48].
Averaging over multiple different measurements of the diffusion at different depths lead to a relatively constant value of 0.047 ± 0.01 μm 2 /s. Importantly this averaged value was comparable to values of the diffusion coefficient for TEWL and Franz cell diffusion measurements found in literature [23,43]. The results suggested that the barrier function in the SC could be located at different depths while still maintaining the same overall barrier properties. Treating the SC with acetone or tape stripping led to an increased diffusion coefficient, as expected, and a much larger SD. The larger SD suggests that the treatments effect on the barrier properties of the SC have large local variation.
The method and model used in this work shows potential for a wide range of applications by offering the ability to quantitatively measure the spatially resolved diffusion of different substances in a variety of biological and other sample types.

| Human skin
The human skin samples used were obtained from breast reduction and abdominoplasty operations. The experiments performed in this work, involving the use of human samples, were approved by the Regional Research Ethics Committee of Southern Denmark, and adhered to the Declaration of Helsinki Principles (2008). Patient consent for experiments was not required because Danish regulations consider human tissue left over from surgery as discarded material. Information on sample preparation and labeling can be found in the supplementary materials. The skin samples used were all abdominal and were donated from three female donors aged 57, 40 and 38 years. Tape-stripping was performed by placing a piece of adhesive tape (14 mm Scotch book Tape 835, M3) firmly on to the skin sample. Subsequently, the tape was removed in a single continual motion. This was repeated five times per sample. Acetone treatment of the samples was performed by covering the skin surface with acetone and left for 2 minutes before wiping the sample with Kimtech wipes.

| D 2 O addition
A sample holder was made by having a needle 0.50 Â 16 mm (B. Braun, Melsungen) bent 90 and glued to the inside of a μ-dish 35 mm, high glass bottom (Ibidi GmbH, Germany). The syringe within the modified dish could then be connected with a silicone hose (diameter = 2 mm and length = 65 cm) to a syringe placed in a syringe pump (Buch & Holm, Denmark). Skin samples of around 1 Â 1 cm in size were thawed and hydrated with milliQ water. The dermis side of the skin samples were then glued to a cover glass with cyanoacrylate adhesive. For some samples one drop of 1 μM fluorescein was applied to SC and left to label for a couple of minutes. One drop of milliQ water was then added to the bottom of the modified μ-Dish in order to minimize refractive index mismatch. A skin sample was placed in the modified μ-Dish on top of the needle with the SC facing downwards. The needle was placed close to the middle of the skin sample. On top of the skin a cover glass was placed and on top of that a small metal weight in order to minimize sample movement. The μ-Dish was then fixated to a holder and placed at the microscope stage. Before CARS imaging, the tip of the needle was found in the microscope and a region of view close to the tip was chosen for the measurements. The syringe pump was set to deliver 1 mL D 2 O with a rate of 12.99 mL/hour.

| CARS microscopy
CARS images of human skin was imaged using a Leica SP8 (Manheim, Germany). The pump laser used was a picosecond pulsed laser at 849 nm (symmetric OD-stretch at 2380 cm À1 ) and the Stokes laser was also a picosecond laser at 1064 nm. Auto-fluorescence was excited by the pump laser (849 nm) and imaged simultaneously.
The system was equipped with two nondescanned PMT detectors in the forward and two in the epidirection. One for CARS and one for auto-fluorescence. Only the epi-direction was used in this experiment. For CARS images a band pass filter (690 ± 25 nm) was placed in front of the CARS detector and for autofluorescent a band pass filter (475 ± 100 nm) was placed in front of the auto-fluorescence detector. Images were collected using the galvo scanners. A Â40 IRAPO water objective (NA = 1.1) lens was used in all the measurements.
Upon starting the experiment, 10 frames were taken to use as a baseline before the syringe pump was started. The images used for the D 2 O analysis were 512 Â 340 pixels (142 Â 142 μm 2 ) and the scan speed was 1.3 seconds per image.

| Image stabilization
Due to sample movement during addition of D 2 O the acquired images has to be stabilized. This was done using FIJI -ImageJ [49] software and the plugin Image Stabilizer [50]. The intensity at different depths with respect to time was measured using a custom written macro. The resulting CSV files could then be used for further data processing in MatLab.

| Mathematical models and data fitting
The models were fitted to the data with a custom MATLAB script using fminsearchbnd. Graphs and statistical analysis were made using Sigmaplot 12.