2.1.1.

Description of skin models

A four-layer skin model (Fig. 1 ) was generated from an eight-layer model used by Zeng ³⁷ to mimic the anatomical structure and optical properties of skin. The eight-layer model was simplified into a four-layer model comprising one layer above the reticular dermis, a reticular dermis layer, a layer below the reticular dermis, and a subcutaneous fat layer. Our interest lies in the reticular dermis $\left(270\mu \mathrm{m}\text{to}1770\mu \mathrm{m}\right)$ , where the microparticle sensors are intended to be introduced and retained in the tissue much like cosmetic dermal-filling agents or tattoos.⁴⁰ We assumed that all of the skin layers above or below the reticular dermis were homogeneous. The four-layer model was designed as a disk of $4\text{-}\mathrm{mm}$ height and $30\text{-}\mathrm{mm}$ diameter; this radial (lateral) extent of each layer was sufficient to be considered infinitely large compared to the size of the excitation light source and sensor models. The wavelength-dependent absorption and scattering coefficients, scattering anisotropy, and refractive index were extracted from multiple sources based on a review of the literature on tissue optical properties.^{37, 41, 42} The optical parameters of the four-layer model were calculated as a weighted average by adding the products of the optical parameters of the eight-layer model and the corresponding volume factors according to their volume contributions to each layer of the four-layer model.^{37, 43} For example, the calculation for the optical parameters of deep dermis in the four-layer model can be expressed as in Eq. 1:

Fig. 1

Eight-layer model and four-layer model. Stratum corneum (SC, $10\mu \mathrm{m}$ ), epidermis (ED, $80\mu \mathrm{m}$ ), papillary dermis (PD, $100\mu \mathrm{m}$ ), reticular dermis (RD, $1500\mu \mathrm{m}$ ), and dermis $(160\mu \mathrm{m}$ ) all have the same optical properties in the model. Upper blood plexus (UBP, $80\mu \mathrm{m}$ ); lower blood plexus (LBP, $70\mu \mathrm{m}$ ); subcutaneous fat $\left(2\mathrm{mm}\right)$ had only refractive index (1.46) in these two models. $\mathrm{SC}+\mathrm{EP}+\mathrm{PD}+\mathrm{UBP}=\text{upper}$ dermis $\left(270\mu \mathrm{m}\right)$ ; $\mathrm{RD}=\text{dermis}$ $\left(1500\mu \mathrm{m}\right)$ ; $\mathrm{LBP}+\text{dermis}=\text{deep}$ dermis $\left(230\mu \mathrm{m}\right)$ .

Table 1

The optical properties for each layer of the four-layer skin model at each of the seven wavelengths used in simulation.

2.1.2.

Description of microsphere sensor models

The sensors were assumed to follow our previous work on microparticle sensor engineering, where the spheres contain platinum octaethylporphine (PtOEP) and rhodamine isothiocyanate (RITC) as indicator and reference dyes, respectively.^{19, 22} Briefly, glucose oxidase and PtOEP were immobilized inside of porous microparticles, which were then coated with diffusion-limited nanofilms tagged with RITC. Using a green LED for excitation and a fiber-optic spectrometer with a bifurcated fiber bundle, the emission of microparticles attached to a microscope slide within a flow-through channel was measured during exposure to phosphate buffers containing different glucose concentrations (flow rate: $4\mathrm{mL}∕\min$ ; oxygen concentration: $277\mu \mathrm{M}$ ). Spectra collected at 0, 50, 100, 150, and $250\mathrm{mg}∕\mathrm{dL}$ were recorded for use in this theoretical study.

The luminescence properties of the microsphere sensors are shown in Fig. 2 . The absorption coefficient $\left(11.16{\mathrm{mm}}^{-1}\right)$ was calculated from the molar extinction coefficient of PtOEP at $540\mathrm{nm}$ , assuming a uniform concentration of $10\mathrm{nM}$ inside the particles (determined experimentally), while scattering internal to the spheres was neglected.²⁰ A refractive index of 1.5 was assumed for the particles, based on their silicate nature. As the size of the particles is an experimentally controllable variable that affects sensor response as well as host response, this parameter was also varied in the simulations to reveal the effects of particle size on emission properties.

Fig. 2

Excitation and emission spectra of sensor particles, as well as spectral properties of a model light source (green filter). For the absorption spectrum shown, [PtOEP] is $10\mathrm{nM}$ . Blue circles and numbers indicate the wavelengths selected for the simulation. (Color online only.)

In this study, the implanted sensors are modeled as “patches” comprising uniform microspheres positioned within the “dermis” of the four-layer model; the microspheres were assumed to be hexagonally packed into a cylindrical slab. Since the particles can be engineered to have different size, and the in vivo packing density can be varied to a degree by controlling the injection concentration, the diameters of microspheres and center-to-center spacing for adjacent spheres (Table 2 ) were varied from $20\mu \mathrm{m}\text{to}100\mu \mathrm{m}$ to observe the influence of these parameters on escaping emission signal. We note that, in reality, one has very little control over the final particle distribution in tissue; however, our studies were aimed at identifying the potential negative consequences of uneven or inconsistent distribution on collection of emission.

Table 2

Geometric parameters of the six microsphere sensor models used in simulations. The sensor names in the simulations are D100S100, D60S100, D20S100, D60S60, D20S60, and D20S20, where the values following the D and S correspond to the diameter and the center-to-center spacing between two adjacent spheres for each case given below, respectively.

2.1.3.

Description of Monte Carlo 3-D modeling system

An example of the entire total simulation system is illustrated in Fig. 3 , including the geometric layout and the excitation source at normal incidence relative to the skin surface. The CAD models the sensor spheres embedded within the skin model; thus, the optical properties of the tissue surrounding and between the spheres were the same as those of the skin model. Each sphere was modeled as an independent object embedded in the skin model, and the set of spheres was arranged into a cylindrical array with different packing efficiency.

Fig. 3

Scheme of simulation system, indicating relative position of excitation light source, implanted luminescent sensing microparticles, and photon capture film. The long-pass filter below the photon capture film is to block scattered excitation light. The distance between the photon capture film and the surface of skin model is $10\mu \mathrm{m}$ , and that between the long-pass filter and the surface of skin model is $9.8\mu \mathrm{m}$ . Although not shown here, the lateral extent of the skin model and the photon capture film $\left(30\mathrm{mm}\right)$ was much larger than the sensor patch size in all cases $\left(1.2\text{to}6\mathrm{mm}\right)$ .

To simulate luminescence emission escaping from the surface of the skin, one excitation wavelength $\left(540\mathrm{nm}\right)$ and six discrete emission wavelengths (570, 585, 620, 635, 645, and $665\mathrm{nm}$ ) were selected from the excitation and emission spectra of sensors to capture the key features (Fig. 3). In the simulations, excitation light undergoes absorption and scattering while traveling through the superficial tissue before striking a sphere. Once the excitation photons enter a sphere, luminescence occurs with isotropic directionality. The luminescence inside a sphere is generated along each excitation ray according to a randomized distance, sampled from a uniform distribution with an average step size between luminescence events that was specified as the radius of the sphere. The total luminescent illumination at each step for each excitation and each emission wavelength is calculated using Eq. 2:

The luminescent light emitted by a sphere, once it escapes the sphere, will undergo further absorption and scattering in the skin model. When rays are incident upon an interface between tissue and surface of another sphere, reflection occurs. Since the absorption coefficient of spheres is high, the luminescence emission from a sphere tends to emit more from the surface than that from the inner space. Moreover, since the excitation light comes from the top of the skin model, the emission tends to be more from the sphere surfaces closest to the skin surface. The luminescent emission rays propagating through the skin model and finally escaping from the surface of the skin model were collected by a photon capture film detector, which recorded the energy and location information of every escaping ray hitting the film. Below the photon capture film, an ideal long-pass filter was attached to remove all excitation light from the captured rays. The collimated excitation light source was perpendicular to the surface of the skin model and was placed directly below the long-pass filter. The photon capture film and the long-pass filter were placed above the surface of the skin model. The distance between the photon capture film and the surface of the skin model was small enough to ensure that the film would capture all the escaping rays. Thus, the photon distribution on the film is an accurate representation of the predicted escaping photon distribution on the surface of the skin model.

Although our microsphere models were not really homogeneous, they were uniformly packed; moreover, the diffusion area of the unit excitation light in the dermis layer was much larger than the cross-sectional area of a unit sphere. Therefore, the output from a narrow incident beam was considered as the impulse response, and convolution was performed to obtain the output response from larger light sources. This approach was validated by comparison of the convolution results and the direct simulation results. The output intensity–spatial distribution of convolution (unit $\text{light}=\mathrm{10,000}$ rays) results and simulation $(\mathrm{12,738}\text{rays}∕{\mathrm{mm}}^2)$ results have no significant difference ( $p\text{-value}=0.9094$ , paired $t$ -test). In OptiCAD, the light source size could not be defined as a pencil beam with infinitely small size, so we used 2-D convolution rather than 1-D convolution. The unit light source was defined as a $100\mu \mathrm{m}\times 100\mu \mathrm{m}$ square, and then the corresponding simulation output was convolved with a larger light source profile. For this work, the light source was considered to be a flat circular incident beam (top-hat) centered at the origin. It contained 10,000 to 4 million rays with varying light source sizes to achieve acceptable variances.

2.2.1.

Assessment of excitation fluence

In some simulations, we assumed that the skin model and the sensor model are infinitely large compared to the unit light source ( $100\mu \mathrm{m}\times 100\mu \mathrm{m}$ square). To justify this, as well as to determine a reasonable dimension for each 3-D model, the depth-dependent spatial distribution of excitation fluence in the skin model was investigated. The problem was broken into two parts: (1) determining the irradiance distribution of a unit excitation light for different depths in tissue; and (2) determining the circular areas that confine different fractions of unit excitation radiant power at different depths in skin model. To extract this information, the skin model was broken up into a series of sections with increasing thickness, and a photon capture film was positioned at the bottom of every skin model section to capture all the exiting excitation light after undergoing absorption and scattering. (The backscattering light reaching each depth from the deeper depth was ignored in these models.) In this situation, the refractive index of the skin model and that of ambient were matched to eliminate the total internal reflection. Thirty simulations were performed in triplicate for each of 30 different depths in the skin model, starting at $30\mu \mathrm{m}$ (close to the upper boundary of upper dermis) and moving to $2\mathrm{mm}$ (close to the lower boundary of deep dermis). In each simulation run, 90,000 rays were used.

2.2.2.

Effects of varying implantation depths

To implement in vivo measurement, it is possible that the implantation depth will be difficult to control with high precision; it is likely that particles will be distributed over different depths in the tissue. Thus, an assessment of measurability versus implantation depth is essential to understand the potential quantitative impact of this practical issue on emission. For these simulations, a single layer of tight-packed sensors $100\mu \mathrm{m}$ in diameter with $100\mu \mathrm{m}$ center-to-center spacing (D100S100) was used. According to the assessment of spatial distribution of excitation fluence in the skin model, the sensor patch was set to be $6\mathrm{mm}$ in diameter (3260 spheres) to ensure that the sensor patch size is much larger than that of the unit input light source ( $100\mu \mathrm{m}\times 100\mu \mathrm{m}$ square). This was verified, as described in Sec. 3.1. The absorption coefficient at $540\mathrm{nm}$ of microparticle sensors was kept at $11.16{\mathrm{mm}}^{-1}$ , assuming a uniform concentration of PtOEP $\left(10\mathrm{nM}\right)$ . The depth of the sensor patch was varied from $400\text{to}1700\mu \mathrm{m}$ , values representative of typical values for the upper and lower boundaries of the dermis layer. Fourteen simulations were performed for each of 14 depths over this range, and each simulation was repeated three times.

2.2.3.

Spatial distribution of escaping luminescence

To design a highly efficient matched optical system for excitation delivery and emission collection, it is necessary to investigate the spatial distribution of the output signal. Using the simulation results for depth-dependent simulations with the $100\mu \mathrm{m}\times 100\mu \mathrm{m}$ unit light source in the previous section, 2-D convolutions were performed to estimate responses to a large flat circular light source with the diameter ranging from $1\mathrm{mm}\text{to}10\mathrm{mm}$ . Convolution was performed with MATLAB.

2.2.4.

Investigation of absolute luminescence response versus different sensor models

To maximize the luminescence output and minimize the number of microparticles and amount of dye for in vivo implantation, it is critical to assess the luminescent output for different sensor configurations and dye distributions. This task can be broken into three steps:

Table 3

Overview of simulations used to investigate sensor size/spacing/layer/effects.

For step 1, the impact of sphere packing efficiency was investigated by varying center-to-center spacing with a constant sphere size. Simulations were completed for three different one-layer sensor models with a sphere size of $20\text{-}\mu \mathrm{m}$ diameter, where center-to-center spacing was varied from $20\mu \mathrm{m}\text{to}100\mu \mathrm{m}$ (D20S100, D20S60, and D20S20). In these simulations, the absorption coefficient of microparticle sensor models was fixed at $11.16{\mathrm{mm}}^{-1}$ ; thus, the effective dye concentration [the product of dye concentration per sphere and sphere packing efficiency (Table 2)] decreased with decreasing sphere packing efficiency (increasing spacing). The maximum diameter of the sensor patches was $1.2\mathrm{mm}$ , and the light source was a $2.4\text{-}\mathrm{mm}$ -diam uniform flat circular incident beam, ensuring that all the sensor particles could be excited by excitation light.

For step 2, the impact of dye concentration was assessed with all six sensor models (D100S100, D60S100, D20S100, D60S60, D20S60, and D20S20). The size of sensor patch models and the light source size were the same as those used in step 1. The absorption coefficients of microparticle sensor models were scaled according to the volume of each sphere and particle numbers so as to make the product of the absorption coefficients, the volume per sphere, and the particle number constant. Thus, each sensor model had the same quantity of dye, and the sensor models with the same center-to-center spacing have the same effective dye concentration.

For step 3, the effects of having more particles in a stacked (layered) structure were assessed with three different multilayer sensor models (D100S100, D60S100, and D20S100). In these simulations, all of the sensor patch models were $4\mathrm{mm}$ in diameter, and a unit input light source ( $100\mu \mathrm{m}\times 100\mu \mathrm{m}$ square) was used. The number of stacked layers of sensors was varied from one to three, and the absorption coefficients of microparticle sensor models were again scaled according to the volume per sphere. Thus, the product of absorption coefficients and sphere volume was constant, and the effective dye concentration in the tissue was constant for different sensor models.

2.2.5.

Simulation for sensor performance after implantation

To investigate the potential to accurately measure changes in sensor implant spectral properties, we investigated the impact of the absorption and scattering of tissue, as well as the microparticle size distribution and concentration for different sensor models with different environmental glucose concentrations. Unless otherwise specified, the absorption coefficient was kept at $11.16{\mathrm{mm}}^{-1}$ (at $540\mathrm{nm}$ ) for all models. We used emission spectra recorded from in vitro analysis of sensors at environmental glucose concentrations of 0, 100, and $250\mathrm{mg}∕\mathrm{mL}$ . Simulations were performed as before, with the six output wavelengths now weighted by the relative emission for different glucose concentrations. The ratio of emission power of $645\mathrm{nm}$ to that of $585\mathrm{nm}$ was calculated from the output spectrum of each sensor emission spectrum corresponding to each glucose concentration using Eq. 3:

The percent change in $645\mathrm{nm}∕585\mathrm{nm}$ emission peak ratios versus glucose concentration was calculated using Eq. 4:

The percent changes of peak ratios were plotted as a sensor response curve. Response curves for six different sensor models were compared. The diameter of the spheres was set to 100, 60, and $20\mu \mathrm{m}$ , and the center-to-center spacing was set to 100, 60, and $20\mu \mathrm{m}$ . The packing configuration of the sensor patch model was maintained at hexagonal packing, and the total number of spheres was constant, as this would best represent the situation of injecting a fixed number of particles with different final distributions. All the sensor models used in this phase were one-layer versions, fixed at an implantation depth of $700\mu \mathrm{m}$ . The light source was the $100\mu \mathrm{m}\times 100\mu \mathrm{m}$ square. The simulation details are given in Table 4 .

Table 4

Models used in assessing effects of sensor configuration on tracking changes in emission.