Sinusitis is an inflammation of the paranasal sinus mucosa, which is treated with antibiotics and often followed by one week of sick leave. Every year, more than 37 million people in the U.S. are diagnosed with sinusitis. Today, the diagnosis of sinusitis is mostly based on the anamnestic history of the patient, and in some cases on paraclinical investigations such as x-ray, ultrasound, and low-dose computerized tomography. Unfortunately, unnecessary antibiotic treatment is very common due to the difficult diagnoses. Therefore, there is a great need for simple, nonintrusive alternatives or complementary methods to detect sinusitis.1, 2
We describe computer simulations related to a new spectroscopic method for human sinusitis detection, which was demonstrated by Persson, Svanberg, and Svanberg.3 This method is based on the sharp spectroscopic absorptive imprint of molecular oxygen characterizing an air-filled cavity traversed by diffusive light. Narrow-band diode laser radiation is launched through the facial skeleton and traverses the frontal sinus cavities in the forehead or the maxillary sinuses in the cheekbone (see Fig. 1 for sinuses location). The light is then diffusely scattered in the deeper lying tissues with part of the light again traversing the cavities and is again scattered in the facial tissue to reach an external detector. For the maxillary sinuses, an alternative light transmission geometry can be employed, where the laser light is transported through an optical fiber connected to the bucal mucosa, scattered through the sinus cavity, and detected externally.
In Ref. 3 we successfully demonstrated sinus cavity monitoring on a healthy human volunteer. The study was supported by extensive experimental trials on tissue-like plastic phantom materials separated by an air gap to simulate the human anatomy.
While all the results could phenomenologically be well understood, we have now performed detailed simulations of light propagation to further guide the development of the technique. For this study, we used the Advanced Systems Analysis Program software, which operates on the Monte Carlo concept. We present simulation data in excellent numerical agreement with experimental data, as well as guidance for improved point monitoring and imaging applications of the technique.
In the next section, we describe the tools used in the study, which include the experimental setup, the computer simulations, and the sinus models investigated in the latter. Then we compare simulations with experimentally obtained data, and we investigate the influence of the scattering coefficient and the detection aperture and the possibility of imaging of the frontal sinuses using the present technique. Finally, conclusions are drawn.
Tools Employed in the Study
The experimental setup used in the measurements and previously reported3 is based on a diode laser operating around where the molecular oxygen A-band is located. The laser was tuned across a sharp individual line [typically a half width at half maximum (HWHM) of a few gigahertz] by applying a saw-tooth ramp to the driving current of the diode laser. To achieve sensitive lock-in detection, a sine-wave was superimposed to the saw-tooth ramp in a wavelength-modulation scheme.
The light was guided to the sample by a fiber. Two detection geometries were used, as shown in Fig. 2 . In the backscattering geometry, a right-angled prism located in the center of the detector was used to launch the light into the scattering medium. An annular aperture with a central block was used to collect photons backscattered from the sample. In the transmission geometry, the fiber was placed on the opposite side of the sample with respect to the detector. The detection was made by a photomultiplier to provide efficient photon collection. A blocking colored glass filter was placed in front of the detector to reduce ambient light.
The signal from the detector was split into two parts: one directly connected to an oscilloscope, and an other processed in a lock-in amplifier. The peak-to-peak value of the oxygen absorptive signal, recorded as the second derivative, was estimated and normalized to the amount of light reaching the detector. By applying a method called standard addition, a so-called equivalent mean path length, , could be estimated corresponding to the distance light has to travel in ambient air to obtain the same signal. This is described in Eq. 1, where corresponds to the gas concentration in ambient air, to the mean path length in the sample, and to the gas concentration in the sample.
The Advanced Systems Analysis Program, , is a ray-tracing software that can be used to simulate light propagation in optical systems. Its implementation is based on the Monte Carlo method and the Henyey-Greenstein approach, which makes it well suitable for simulations in biological tissue and other light scattering materials.4 The basic idea of Monte Carlo is to apply a probability view of light propagation, where an individual photon can be absorbed or laterally scattered by a certain angle according to relevant probability laws applying to the situation. After launching a large number of individual photons, the probability distributions are transformed into intensity distributions that can be compared to measured data. Basically, the physical experiment is repeated in the computer, considering a multitude of identical incident photons that later experience different fates.5
All simulated models considered in the present work incorporate scattering materials representing the human tissue and an air space representing the sinuses. By studying the rays reaching the detector in the simulations, the equivalent mean path length could be calculated for all models. For every ray, computes information such as the flux at every surface it enters or is reflected off, as well as the distance it has traveled in all media. We then calculate the equivalent mean path length by adding up, for all rays, the product of the distance each individual ray has traveled in the air space and its corresponding flux when reaching the detector, and dividing it by the total flux reaching the detector. This is described in Eq. 2, where is the total numbers of rays reaching the detector.
It should be noted that rays reaching the detector that have not traveled in the air gap will have , and only contribute to the simulated equivalent mean path length with a diluting factor. These rays are referred to as short-cut photons.
Investigated Sinus Models
The phantom model representing the human sinuses used in the experiments and the simulations contains two scattering slabs made of white Delrin-type plastic, with optical properties similar to human tissue (see Fig. 2). The scatterer closest to the light source is referred to as the primary scatterer , and the other one is referred to as the secondary scatterer . In the case of the transmission geometry, is the scatterer closest to the detector, and in the backscattering geometry it is the one located above . The thicknesses of the scatterers are denoted and . In-between the scatterers, a variable air gap was placed. The annular aperture determining the detection area was simulated by a variable absorptive mask.
To describe imaging of the frontal sinuses, a new model was created for simulations. In the case of imaging of the frontal sinuses, the sinuses were represented by two ellipsoids embedded in a scattering medium. The same backscattering detection geometry shown in Fig. 2a was used.
Results and Discussion
Comparison between Experiments and Simulations
To simulate the phantom experiments described in Ref. 3, which represent measurements on the human sinuses, we created a model in consisting of two scattering materials and , respectively, with a variable air gap in-between (see Fig. 2). The thicknesses and were varied, consistent with the experiment. The optical properties (the anisotropy factor , the scattering coefficient , and the absorption coefficient ) were measured for one of the Delrin plates used in the experiments by the integrating-sphere technique6, 7 to be , , and at . These values were used in the simulations. A small absorption was added to better fit the experimental data. The index of refraction for the Delrin plates was set to 1.48 and for the air gap to .
Initial simulations, including all reflection and refraction phenomena according to the Fresnel and Snell formulas, were performed. It was noticed that the influence of the reflection was of minor importance for angles below the critical angle in our case. Thus it was suppressed in subsequent simulations, resulting in a manageable computer time. Above the critical angle, the total internal reflection was fully taken into account. The motivation to this approach has its roots in the fact that the solid angle for impinging photons corresponding to total internal reflection is much larger than that of shared reflection and refraction. Further, for photons impinging toward an optically denser medium, reflection has a small contribution as compared to penetration into the medium for the dominant fraction of the solid angle. For example, with two 10-mm-thick Delrin plates spaced by , the underestimation of the true is about 8% and will, for larger spacings, become smaller (numerical values are provided in the Appendix.
In Fig. 3 the simulated and the experimental are shown for the backscattering detection geometry, representing measurements on the frontal sinuses [model according to Fig. 2a]. The inner detection radius was about , and the outer about in both simulations and experiments. Each simulated data point consists of launching rays into the model.
As during the experimental trial,3 it was observed in the simulation that an thicker than did not influence the outcome. This can be understood because the photons being scattered deeper than into have a very low probability to be detected. As can be seen in Fig. 3, the agreement between the simulations and the experiments is very good. It can be seen that introducing a thicker secondary scatterer results in increasing. This is because a thicker secondary scatterer results in a higher probability for the photons to be scattered back to the primary scatterer and be detected after passing through the air gap. On the contrary, a thicker primary scatterer results in a decrease of . The photons have now a lower probability of reaching the secondary scatterer, and even if they do, the backscattered photons from the secondary scatterer, are now more likely to be scattered outside the finite size of the detector. Because of the finite size of the detector, the signals first increase with increasing air gap until a certain distance, where they start to decrease. The maximum of the signal should be located at an air gap distance as large as possible to prevent two air gaps from giving rise to the same value, which could be critical in the case of real human measurements. The signal behavior is explained in more detail in Ref. 3.
The deviation between the experimental data and the simulated data can be due to many reasons. In the experiments, the thicknesses of the Delrin plates, the air gap, and the annular detection aperture could not be determined as precisely as they can be in the simulations. During the simulation it was also observed that a small change in one of the optical properties influenced the outcome. It should also be noted that the same optical properties were used for all the different Delrin plates in the simulations. This might not be the case, since the plates used in the experiments did not come from the same batch.
Experimentally it is impossible to explore the influence of multiple passes through the air gap. In contrast, this can easily be done in the simulations. Table 1 shows the number of times the detected photons pass over the air gap, for and , and for different air-gap distances. It should be noted that only even numbers of passes are possible, since the photons have to travel back and forth to be detected. Similar results were observed for the other thicknesses discussed in Fig. 3.8
Distribution of the number of passes through the airgap for l1=10mm and l2=30mm and different air-gap distances in the backscattering detection geometry according to Fig. 2a.
|Air distance [mm]||1||5||10||30||70|
|2 passes [%]||31||42||51||74||98|
|4 passes [%]||20||22||23||19||2|
|6 passes [%]||14||13||11||5||0|
|8 passes [%]||10||8||6||1||0|
|10 passes [%]||7||4||3||1||0|
|12 passes [%]||5||4||2||0||0|
|14 passes [%]||3||2||2||0||0|
It can be seen that signals due to multiple passes can be significant, in contrast to what was intuitively assumed in Ref. 3. This illustrates the power of detailed Monte Carlo simulations, taking into account the actual regime of values of the optical constants.
It can also be seen from Table 1 that a larger air-gap distance results in fewer multiple passes. This phenomenon can be explained by the finite size of the detector. It was observed from the simulations that a thicker resulted in a increase of multiple passes. This can be understood because a thicker primary scatterer can collect oblique rays to be detected by the finite size detector that have been through more multiple passes. It was also observed that a thicker secondary scatterer increases the number of multiple passes, as a thicker secondary scatterer allows more photons to be scattered back to the primary scatterer.8
In Fig. 4 the simulated and the experimental are shown for the transmission detection geometry, representing measurements on the maxillary sinuses [model according to Fig. 2b]. The mask radius was about in both simulations and experiments. Each simulated data point consists of launching rays into the model.
Again, very good agreement between the experimental and simulated data can be seen. The obtained for the different scatterer parameters increases with increasing air-gap distance. This can be understood, since in the case of the transmission detection geometry, all photons detected must have traveled over the air gap. This is in contrast to the backscattering detection geometry, where the short-cut photons (photons that have only been scattered in before being detected) dilute the observed . The signal behavior is explained in more detail in Ref. 3.
It can be seen that the photon statistics are not as good as for the backscattering detection geometry. For the transmission geometry, the percentage of detected photons in relation to the number of incident photons is much lower than for the backscattering case, despite the fact that all photons are carrying information of the sinuses. A much longer simulation time is needed to achieve good photon statistics. Again, the deviation between the experimental and the simulated data can be due to many reasons, as discussed before.
The number of multiple passes was also investigated in the transmission detection geometry for the models representing measurements on the human maxillary sinuses. Table 2 shows the number of times the detected photons pass the air gap, for and , and for different air distance. In this geometry, only odd number of passes are possible. Similar results were observed for the other thicknesses discussed in Fig. 4.8
Distribution of the number of passes through the air gap for l1=10mm and l2=10mm and different air-gap distances in the transmission detection geometry according to Fig. 2b.
|Air distance [mm]||1||5||10||30||70|
|1 passes [%]||35||49||61||81||88|
|3 passes [%]||22||23||20||14||10|
|5 passes [%]||14||11||9||4||2|
|7 passes [%]||9||6||4||1||1|
|9 passes [%]||6||4||2||0||0|
|11 passes [%]||4||3||2||0||0|
|13 passes [%]||3||1||1||0||0|
Again, as in the case of backscattering detection geometry, it can be seen that signals due to multiple passes can be significant. It can also be seen that for a larger air-gap distance, fewer multiple passes occur. As discussed before, this phenomenon can be understood from the finite size of the detector. It was also observed that a thicker results in more detected multiple pass rays. This can be understood, since more oblique rays will be collected by the finite size of the detector. The inverse argument is applicable to explain why a thicker results in less multiple passes.8
Influence of the Scattering Coefficient
The scattering coefficient and the thickness of the facial bone may vary between different humans, especially with age.9 The scattering coefficient of bone is significantly smaller for children than for adults. In our model, a lower scattering coefficient has similar influence as a smaller thickness of the primary scatterer, i.e., a lower scattering coefficient of results in a higher probability for the photons to reach be scattered back and be detected by the finite size of the detector.
Figure 5 shows the influence of the scattering coefficient of and . The simulations were performed for the backscattering detection geometry shown in Fig. 2a, with an annular detection aperture of dimensions and . The optical properties used in the simulations are taken from different sources and presented in Table 3 .10, 11, 12 In all simulations, except for Delrin, was used. For each trial with a set of optical properties, simulations were performed with the same air-gap intervals as in Figs. 3 and 4, and by launching rays into the model for each air-gap distance.
Optical scattering properties used in the simulations taking from different sources.10, 11, 12
|Optical properties||μs [mm−1]||μs′ [mm−1]||μa [mm−1]||g||λ [nm]|
|Adult skull in vitro||16||2.0||NA||0.87||800|
|Adult skull in vivo||13||0.9||NA||0.93||674–849|
|Neonate skullin vivo||9.3||0.65||NA||0.93||761|
From the simulations it was observed that only the reduced scattering coefficients and influence the result. The difference between optical properties corresponding to adult skull in vivo and neonate skull in vivo is only the reduced scattering coefficient. By comparing the curves obtained for these cases in Fig. 5, it can be seen that a lower results in a larger air-gap distance needed to obtain the maximum signal. This is because a decrease in can also be seen as a decrease in the thickness of the primary scatterer. This phenomenon can also be seen in Fig. 3a. The influence of the absorption coefficient can be seen by comparing the simulated curves obtained for Delrin and adult skull in vitro. A lower gives rise to a higher integrated signal.
Annular Detection Aperture
In the experiments, the choice of annular detection aperture was done intuitively. Here we investigate the influence of the aperture geometry using the simulations.
In Fig. 6 , is plotted as a function of the air-gap distance for different annular detection apertures in the case of the backscattering detection geometry. The thicknesses of and were set to 5 and , respectively. The simulations were done with neonatal optical scattering properties (shown in Table 3). Again, the simulations were performed with the same air-gap intervals as in Figs. 3 and 4, and by launching rays into the model for each air-gap distance.
In Fig. 6a, it can be seen that increasing while keeping constant increases the signal. This is because a larger allows more photons traveling over the air gap to reach the detector. This also gives rise to a larger photon flux reaching the detector. In Fig. 6b, is kept constant while is varied. It can clearly be seen that increasing results in a higher . This is equivalent to blocking out more of the short-cut photons that have only been scattered in the first medium, and thus only influence by diluting it. As can be seen in the figure, an increase in also pushes the maximum of the signal toward larger air-gap distances. This might be very important for preventing two different air-gap distances giving rise to the same value.
The conclusion from these different curves obtained with different annular detection apertures is that an outer detection aperture as large as possible should be used. In the current setup, a photomultiplier (Hamamatsu 5070A) with a radius of is used, which limits the outer radius. The simulations also show that an inner radius as large as possible should be used with an upper limit set by the experimental signal-to-noise ratio. Increasing the power of the diode laser used should thus be advantageous. Similar results were obtained and the same conclusion was drawn for scattering properties according to an adult, and are discussed in detail in Ref. 8.
A study of the influence of the detection aperture in the transmission detection geometry was also performed. In this detection geometry, the size of the aperture did not have the same influence on the signal as in the backscattering detection geometry. As mentioned before, all photons reaching the detector have traveled over the air gap distance, and short-cut photons need not to be blocked. This is discussed in more detail in Ref. 8.
To investigate the feasibility of imaging the human frontal sinuses, a model according to Fig. 7a was studied. The two frontal sinuses were represented by two air-filled ellipsoids ( high, wide, deep) separated by and embedded at a depth of in a scattering material, with optical properties corresponding to adult in vitro (Table 3). The backscattering detection geometry used in the simulations is shown in Fig. 2a. Figure 7b shows the results of the simulations when using an annular aperture with and . The whole image consists of 1000 simulated data points with the light and detection arrangement scanned across the model. For each position, rays were launched into the sample. The locations of the ellipsoids are also indicated. The figure clearly shows the possibility of spatially resolving the cavities.
Many factors might be considered to determine what constitutes a “good” image. It might be important to maintain the largest possible contrast between maximum and minimum signals, e.g., corresponding to the centers of the ellipsoids (the sinuses), and the center gap between them. To be able to resolve the two sinuses, it might also be critical to use a sufficiently small detection aperture to not blur the image.
In Fig. 8 , the cross sections of the simulated backscattered image are shown for different detection apertures. In Fig. 8a, was kept constant at , while was varied between 5 and . The real cross section is also included for the sake of clarity. It can be seen that selecting a larger outer detection radius results in a higher value. In Fig. 8b, was kept constant at , while was varied between 4 and . It can be seen that the spatial resolution is rather independent of the chosen detection radii. Then it is important to obtain a sufficiently large value while maintaining a large photon flux.
Numerical values for the phantom simulations when employing the backscattering detection geometry (l1=l2=d=10mm) . The different cases are explained in the text. For each case, 0.5∙106 rays were launched into the sample. In case 3, the criterion that each ray was allowed no more than two divisions was set.
|Case 1||Case 2||Case 3|
|Total number of detected rays||8670||59163||67676|
|Detected oxygen rays||22||1788||3188|
Conclusion and Outlook
The power of Monte Carlo simulations in understanding and optimizing the experimental conditions for human sinus monitoring is demonstrated. The possibility of working with arbitrary geometries in the approach is particularly valuable. The simulations agree well with experimental data for the limited setup conditions studied so far. This gives us confidence in using the simulations presented in the process of developing a new and clinically relevant modality for human sinus diagnostics.
In particular, realistic insight in the imaging capability for the frontal sinuses is obtained, showing that the main features of the cavities can be retrieved noninvasively. The geometry is complex and more difficult to test experimentally, but can be illuminated in the simulations. Since we deal with highly scattering media, only an equivalent air path is retrievable from the experiments, but this parameter seems to relate well to the factual geometry. Insight into how sensitive the oxygen signals are to the detailed probe positioning is obtained, which is valuable in nonimaging applications. Among these, the measurement of gas-dynamic changes through the channels connecting the sinuses to the nasal cavity, and the measurement of concentration ratios between two gases, such as oxygen and water vapor, are of particular interest, since they yield information on the conductivity and absolute oxygen concentration, independent of the scattering properties.
The authors would like to thank Mats Andersson for technical assistance, and Breault Research for providing the software. We are also grateful to Katarina Svanberg for help with the medical aspects of this work, and to Stefan Andersson-Engels for valuable comments and advice. This research was supported by the Swedish Research Council, and the Knut and Alice Wallenberg Foundation.
In the simulations presented, Fresnel reflection below the critical angle was not taken into account. It was observed that the influence was of minor impact in initial simulations, and therefore time consuming in the simulations and could be avoided. For example, 200 injected rays result in about 200 000 rays to track when simulating the backscattering detection geometry for the case of 10-mm-thick Delrin plates spaced by , including Fresnel reflection. This is due to the fact that at each surface, the rays get divided up into one transmitted ray and one reflected ray with different weights according to Fresnel’s laws. In Table 4 , numerical values are given for the phantom simulations when employing the backscattering detection geometry. Data are presented for the case of 10-mm-thick Delrin plates (index of reflection equal to 1.48) spaced by (index of reflection equal to 1). In the table, the total number of detected rays, the number of rays detected that have traveled over the air gap (oxygen rays), and the resulting simulated values are shown. Numerical values from three different cases are presented.
Case 1. No Fresnel reflection is taken into account in the simulation, but refraction according to Snell’s law is included.
Case 2. Total internal reflection is taken into account in the simulations, but Fresnel reflection is excluded below the critical angle. Snell’s law refraction is included.
Case 3. All reflection and refraction phenomena according to the Fresnel and Snell formulas are taken into account in the simulations.
It can clearly be seen that by not taking Fresnel reflection into account (Case 1), no detectable simulated is obtained. It can also clearly be seen that the important effect, which needs to be taken into account in the simulations, is total internal reflection (Case 2). By comparing Case 2 with a “correct” (and extremely time-consuming) calculation (Case 3), an error on the of only about 8% when excluding Fresnel reflection below the critical angle from the simulations is obtained.