How to model the human skin and optical heart rate sensors in OpticStudio

Also available in 日本語, 简体中文

Photoplethysmography (PPG) is a low cost, non-invasive optical technology that takes physiological measurements on the surface of the skin. One of its most widespread applications is the wearable heart rate sensor included in commercially available smartwatches and sport bracelets, that provides comfortable and continuous pulse monitoring during everyday tasks. This article demonstrates how to model the human skin in Zemax OpticStudio for physiological measurements, and illustrates a time-dependent simulation of PPG-based heart rate sensors using ZOS-API.

Authored By Csilla Timar-Fulep

download example (.zip)

Introduction

PPG devices consist of infrared or visible range light-emitting diodes (LEDs) and photodetectors. They provide a simple optical technique to detect blood volume changes in tissues. Light is more strongly absorbed and scattered by blood than by the surrounding tissue. Therefore, the pulsation of blood causes a variation of opposite phase in the signal of the detector. This article describes how to implement a human skin tissue model in OpticStudio, and shows how to simulate the measured signal of a PPG device over time using a ZOS-API application.

Basic design

PPG sensors can be used either in reflection or in transmission mode. As the penetration depth of the light depends on its wavelength, green and yellow LEDs are most suitable to take measurements in the superficial blood flow and are typically used in reflection mode. On the other hand, infrared and near-infrared wavelengths are better suited for measurements of deep-tissue blood flow and can be used in transmission mode as well. In this use case, we present a reflection PPG device.

Our goal was to develop a realistic skin model based on data published in the relevant literature. Thus, we intended to apply a wavelength at which the optical parameters of both the skin and the blood are widely available in the literature, and which is also close to the most frequently used wavelengths in commercial devices. As a result, we chose the modelling wavelength of 575 nm, and we used a QSMF-C160 LED (Avago Technologies) as the source, modeled as a Source File in Zemax OpticStudio.

Human skin model

To simulate light transport in tissue media, we implemented a layered skin model, which accounts for the epidermis, the dermis, and for the subcutaneous fat. As the primary goal for this use case was to simulate a PPG-based heart rate sensor, where the key point is to measure the changes caused by the pulsation of the blood, we focused on accurately modelling the layers where this pulsation can be observed. Therefore, we separately modelled the sub-layers of the dermis with different blood content values, namely the papillary dermis, subpapillary dermis, upper blood net dermis, reticular dermis, and deep blood net dermis. On the other hand, as there is no blood content in the epidermis - to keep the model simple - we decided to use only one thick epidermis layer, which accounts for all the stratum corneum, stratum granulosum, stratum spinosum, and stratum basale. Finally, similarly to most of the published skin models, we represented the subcutaneous fat with one single layer too.

All the above-mentioned layers are modelled as Rectangular Volumes in OpticStudio. The thickness values of the layers were based on literature data, and the cross-sectional size was determined in such a way that there is no light leakage at the sides. The subsequent layers were placed by using the previous layer as Reference Object and applying Pickup solve on the Z Position value from the Z Length column of the previous layer. This solution ensured that the layers were located closely one after the other without any gap between them.

Customizing the tissue layers

As this case study relies solely on data published in the literature, we did not carry out any new measurements throughout the research process. Although the model parameters were based on published data, it is important to note that the optical parameters of the human skin may vary significantly across the population. Thus, different parameters might be required for specific subject groups. Therefore, if more accurate data are available for your specific applications, then do not hesitate to customize the model accordingly! 

Representing all the individual blood vessels in the skin would require the addition of several hundreds of objects with complex spatial arrangement, and would reduce the generality of the model, thus this type of modelling is not widespread in the literature. Therefore, we did not apply this technique either. Instead, we took the blood content of the different skin layers into account by calculating the weighted average of the optical parameters of the blood and the surrounding tissue.

Accordingly, we modelled the material of the skin layers using the Model material solve type in OpticStudio, based on the following raw data:

Bulk scattering in tissue

Light scattering by small particles in turbid media, such as in biological tissues, can be accurately described by the Henyey-Greenstein distribution function. The Henyey-Greenstein model has only one free parameter, the anisotropy factor g. The domain of this parameter is the [-1, 1] interval, where g=-1 corresponds to backscattering, g=0 describes isotropic scattering, and g=1 represents forward scattering. The angular distribution of the scattered light is defined as:

In the Non-Sequential Mode of OpticStudio, the Henyey-Greenstein bulk scattering model is available in a DLL (Henyey-Greenstein-bulk.DLL), included in the OpticStudio installation files. Further explanations of the model together with comprehensive analyses of the DLL are presented in the knowledgebase article Using the Henyey-Greenstein distribution to model bulk scattering. Besides, a detailed discussion on the surface and bulk scattering models in OpticStudio can be found in this knowledgebase article: What scattering models are available in OpticStudio?

In the attached multilayer skin model, the scattering parameters of each layer were set based on realistic values presented in the literature. While the input parameters of the Henyey-Greenstein scattering DLL are the Mean Path, the Transmission fraction, and the g anisotropy parameter, in the literature, typically the scattering and absorption coefficients, µ_s and µ_a respectively, are published with the anisotropy factor. Therefore, we used the following equations to compute the input parameters of the model:

Similarly, as before for the refractive index, the scattering and absorption coefficients, as well as the anisotropy factor of the different skin layers were calculated as the weighted average of the corresponding values for the blood and for the rest of the tissue, based on the following raw data corresponding to the applied wavelength of 575 nm:

The multilayered skin model and light transport in the tissue is shown in the 3D layout plots below. To illustrate the individual scattering events, rays are colored by segments on the plots.

To provide numerical results in addition to the illustration, we added three Detector Rectangle objects to the design. They are separated by a thin air gap from the surface of the skin. Two detectors have the same cross-sectional size as the skin layers, one facing towards the source and the other towards the skin model, to measure all the incoming and back-scattered light respectively for references. The third detector is a small one (2 mm by 2 mm), again facing towards the skin, which represents a typical photodetector in PPG devices.

The above-presented design can be found in the Article Attachments (skinModel.zar), and it can be used as an off-the-shelf skin model when the time-dependence of the measurements/simulations is not relevant. On the other hand, the way to model time-dependent effects - e.g., in case of heart rate sensors - is discussed in the following section.

Heart rate sensor simulation

To simulate heart rate monitoring, we used the ZOS-API to mimic pulsating blood flow in the tissue. We modelled the different phases of the cardiac cycle by adjusting the blood content of the skin layers, and then we examined the detected back-scattered light as a function of the time steps. Temporal changes in the blood content of the layers were taken into account by a multiplication factor, assuming the blood volume varies proportionally and simultaneously in each layer. In this example, we used the Python API - connected to OpticStudio via .NET - to modify the model parameters, run ray traces with the finetuned settings, and finally to analyze and plot the results. The Python script of the Standalone Application can be downloaded from the Article Attachments (PPGsimulation.py).

Modifying the tissue parameters from the API

According to the literature, as the heart pumps blood into the vessels in the systole phase the relative blood content of the skin layers can double. We used an empirical function to characterize this pulsation, it is illustrated in the image below for 5 cardiac cycles with 10 steps/cycle.

Based on this, first we calculated the blood content in each layer, and then we updated the refractive index, as well as the mean free path, the transmission, and the g anisotropy factors of the Henyey-Greenstein scattering distribution accordingly. The corresponding Python code reads:

layer = TheNCE.GetObjectAt(layerNum) 
solver = layer.MaterialCell.CreateSolveType(ZOSAPI.Editors.SolveType.MaterialModel)
solver._S_MaterialModel.IndexNd = n
layer.MaterialCell.SetSolveData(solver)
volPhysData = layer.VolumePhysicsData
volPhysData.ModelSettings._S_DLLDefinedScattering.MeanPath = meanPath
volPhysData.ModelSettings._S_DLLDefinedScattering.SetParameterValue(0, transmission) 
volPhysData.ModelSettings._S_DLLDefinedScattering.SetParameterValue(1, g)

Finally, we ran a ray trace for each time step, and pulled the results, i.e. the total number of ray hits and absorbed power, from the small PPG detector. To improve ray tracing speed, we applied simple ray splitting, the details of which are discussed in the knowledgebase article What is Simple Splitting?. Ray tracing and data extraction can be accomplished with the following Python commands:

NSCRayTrace = TheSystem.Tools.OpenNSCRayTrace()
NSCRayTrace.ClearDetectors(0)
NSCRayTrace.SplitNSCRays = True
NSCRayTrace.ScatterNSCRays = True
NSCRayTrace.UsePolarization = True
NSCRayTrace.IgnoreErrors = True
NSCRayTrace.SaveRays = False
NSCRayTrace.Run()
NSCRayTrace.WaitForCompletion()
NSCRayTrace.Close()

hits = -3 # pixel =-3 for total hits
power = 0 # pixel = 0 for total power
Data = 0 
hits_bool_return, total_hits = TheNCE.GetDetectorData(detectorNum, hits, Data, 0)
power_bool_return, total_power = TheNCE.GetDetectorData(detectorNum, power, Data, 0)

Results

Since the absorption and scattering coefficients of the blood are much larger than the rest of the tissue, higher blood content causes lower amount of back-scattered light as well as lower measured power on the detectors, and vice versa. Based on our simulations, using a 1 W source and 10⁵ analysis rays in each time step, the modelled pulsation of the blood content results in a ~10…15% variation in the detector signal, which is shown in the plot below.

These results are in good agreement with experimental observations presented in the literature. In commercial heart rate sensors typically a straightforward signal processing algorithm is applied on this data as part of post-processing. As a first step smoothing is performed on the noisy signal, and then the number of peaks above a certain pre-defined threshold are counted per minute to calculate the heart rate.

References

1. Maeda, et al. Monte Carlo simulation of spectral reflectance using a multilayered skin tissue model. Optical Review (2010)
2. Sinichkin, et al. In vivo fluorescence spectroscopy of the human skin: experiments and models. Journal of Biomedical Optics (1998)
3. Meglinski, Matcher. Quantitative assessment of skin layers absorption and skin reflectance spectra simulation in the visible and near-infrared spectral regions. Physiological Measurement (2002)
4. Doronin, et al. Assessment of the calibration curve for transmittance pulse-oximetry. Laser Physics (2011)
5. Meglinski. Monte Carlo simulation of reflection spectra of random multilayer media strongly scattering and absorbing light. Quantum Electronics (2001)

KA-01986