Project 1

Simulation of Light Tissue Interaction

  • Institute: HOT - Hannoversches Zentrum für optische Technologien, Hannover
  • Principle Investigators: Dr. Merve Wollweber, Prof. Dr. Bernhard Roth
  • Researcher: Dr. Oliver Melchert

In subproject 1 the focus is on various numerical aspects of optoacoustics in the context of computational biophotonics. Our aim is to perform numerical experiments and to devise measurement protocols that facilitate a better understanding of the features of optoacoustic (OA) signals that result from melanin enriched absorbing structures within tissue.

In summary, optoacoustics is a two-part phenomenon consisting of

  • Part I: Optical absorption of laser beams inducing photothermal heating of tissue
  • Part II: Acoustic emission of ultrasound waves due to thermoelastic expansion and stress field relaxation

Whereas the optical absorption is assumed to occur instantaneously, the acoustic propagation of sound waves is a comparatively slow process that occurs on a microsecond timescale. In this regard, note that typical propagation distances are on the order of centimeters and that the propagation of acoustic waves in soft tissue (i.e. elastic solids) occurs with a speed of 1400-1600 meters per second. Hence, in subproject 1 the challenge is to combine the absorption of laser light by tissue and the acoustic propagation as a multitimescale problem. Since both underlying phenomena occur on vastly different timescales we can adopt a divide-and-conquer strategy, i.e. decouple the optical absorption problem from the acoustic propagation problem and solve both subproblems independently. Our research activity is centered around three main topics:

  • Topic 1: The direct OA problem
  • Topic 2: The inverse OA problem
  • Topic 3: Algorithms in computational biophotonics

Topic 1: The direct OA problem

Here, our aim is to understand the principal features of optoacoustic (OA) signals resulting from measurements on melanin enriched absorbing structures within tissue, and to numerically verify OA signals generated in mulit-layered PVA hydrogel phantoms with focus on (nondestructive) OA depth profiling.

One of the milestone achievements in HYMNOS subproject 1 was to collaborate with project MeDiOO at HOT in order to complement optoacoustic signals detected in the lab by means of numerical simulations. The figure below (taken from an early draft of our joint article) illustrates the good agreement of the measured (solid orange curves) and numerically predicted (dashed blue curves) excess pressure curves for (a) a three-layered PVA hydrogel based tissue phantom, (b-c) four-layered tissue phantoms.

In our numerical studies we adressed the prediction of OA signals observed in the acoustic near- and farfield in both, forward and backward detection mode, in PVA-H tissue phantoms, see here and here.

Topic 2: The inverse OA problem

Here, our aim is to solve the optoacoustic (OA) source reconstruction problem in order to invert OA signals to inintial stress profiles and to infer the optical properties of the underlying material. Therefore we consider the OA problem in the paraxial approximation where the source reconstruction is achieved by the inverse solution of a Volterra integral equation.

The figure below illustrates an exemplary solution of the OA source reconstruction problem using a Picard-Lindelöff iteration scheme for synthetic input data (solid blue curve). As evident from the figure, the inverted signal (solid red curve) compares well to the true initial stress profile (solid black curve). The intermediate auxiliary stress profiles (dashed gray lines) feature a pronounced rarefaction dip that shifts outwards upon iteration. The iteration procedure is stopped as soon as the Chebychev-norm between two subsequent auxiliary profiles decreases below a fixed threshold.

In our numerical studies we also analyzed the effects of noise, detector-to-sample distance and finite detector size on the OA signals and the reconstructed initial stress profiles, see here.

Topic 3: Algorithms in computational biophotonics

Here, we study efficient algorithms for the accurate forward and reverse evaluation of the discrete Fourier-Bessel transform. The transform is used as a numerical tool facilitating a polar convolution of two radially symmetric functions. This is relevant for applications in tissue optics and optoacoustics, where a recurrent task is to perform a beam-shape convolution in order to yield the material response to an extended laser beam from its Greens-function response. The latter results from a more complex measurement process modeled in terms of a Monte Carlo approach and is, in the worst case, known on a finite sequence of coordinate values, only. The considered applications require the repeated (hundreds of times) calculation of a forward and reverse Fourier-Bessel transform. Thus, time-efficiency is key and fast numerical procedures are valuable.

The figure below illustrates the laser profile convolution procedure for different irradiation source profiles (ISP). Colums from left to right: Gaussian ISP, Flat-top ISP, and Donut ISP. Rows from top to bottom: increasing scattering anisotropy factor.

We applied the developed convolution procedure to the problems of finite-size beam-shape convolution in polar coordinates and the prediction of photoacoustic transients observed in actual experiments, see here (in the linked article we address the important issue of testing research code by illustrating a sequence of unittest).

Developed Software

  • LightTransportMC - A simple Monte Carlo model of photon transport and fluence rate estimation in semi-infinite homogeneous absorbing and scattering media
  • SONOS - A fast Poisson integral solver for layered homogeneous media
  • INVERT - Optoacoustic inversion via Volterra kernel reconstruction
  • PCPI - A python module for optoacoustic signal generation via Polar Convolution and Poisson Integral evaluation
  • dFBT-FJ - Efficient polar convolution based on the discrete Fourier-Bessel transform for application in computational biophotonics
  • LEPM-1DFD - One-dimensional finite-difference code for piecewise homogeneous, linear elastic and piezoelectric materials
  • in collaboration with TP2: BPPA-2DFD - Algorithms and benchmark instances used for beam propagation in the paraxial approximation via 2D finite-difference time domain methods
  • in collaboration with TP4: FMAS - Forward model for the analytic signal in ultrashort pulse propagation