Detection, numerical simulation and approximate inversion of optoacoustic signals generated in multi-layered PVA hydrogel based tissue phantoms
DDetection, numerical simulation and approximate inversion of optoacousticsignals generated in multi-layered PVA hydrogel based tissue phantoms
E. Blumenr¨other, O. Melchert, M. Wollweber, B. Roth
Hannover Centre for Optical Technologies (HOT), Interdisciplinary Research Centre of the Leibniz Universit¨at Hannover,Nienburger Str. 17, D-30167 Hannover, Germany
Abstract
In this article we characterize optoacoustic signals generated from layered tissue phantoms via short laserpulses by experimental and numerical means. In particular, we consider the case where scattering is ef-fectively negligible and the absorbed energy density follows Beer-Lambert’s law, i.e. is characterized by anexponential decay within the layers and discontinuities at interfaces. We complement experiments on sam-ples with multiple layers, where the material properties are known a priori , with numerical calculations fora pointlike detector, tailored to suit our experimental setup. Experimentally, we characterize the acousticsignal observed by a piezoelectric detector in the acoustic far-field in backward mode and we discuss theimplication of acoustic diffraction on our measurements. We further attempt an inversion of an OA signalin the far-field approximation.
Keywords:
Optoacoustics, PVA hydrogel phantom, approximate signal inversion
PACS:
1. Introduction
Recent progress in the field of optoacoustics(OAs) has been driven by tomography and imag-ing applications in the context of biomedical optics.Motivated by their immediate relevance for medi-cal applications, high resolution scans on living tis-sue proved the potential of the optical absorptionbased measurement technique [1, 2, 3]. Requiringa multitude of detection points around the sourcevolume, OA tomography allows for the reconstruc-tion of highly detailed images, see, e.g., Ref. [4], as-suming a mathematical model that mediates theunderlying diffraction transformation of OA signalsin the forward direction [5, 6, 7, 8].However, for most cases of in vivo measurements,especially on humans, it is not feasible to place ul-trasound detectors in opposition to the illumina-tion source (with the “object” in between), i.e. to
Email addresses: [email protected] (E.Blumenr¨other), [email protected] (O. Melchert) work in forward mode . Instead, it is worked in back-ward mode , where detector and irradiation sourceare positioned on the same side of the sample. Us-ing elaborate setups combining the paths of lightand sound waves it is possible to co-align opticaland acoustic focus within the sample. By scanningover a multitude of detection points it is then possi-ble to produce 3D images with very high resolution,see, e.g., Ref. [9].A conceptually different approach is to performmeasurements by means of a single, unfocusedtransducer. Albeit it is not possible to reconstructOA properties of arbitrary 3D objects with a fixedirradiation source and a single detection point only,useful information of the internal material proper-ties of, say, layered samples can be gained never-theless. In this regard, acoustic near-field measure-ments by means of a transparent optoacoustic de-tector where shown to reproduce the depth profileof absorbed energy density and absorption coeffi-cient without the need of extensive postprocessing[10, 11, 12, 13]. However, requiring close proxim-ity and plane wave symmetry, near-field conditionsare unrealistic considering most measurement sce-
Preprint submitted to Photoacoustics November 11, 2018 a r X i v : . [ phy s i c s . b i o - ph ] M a y arios. In contrast, the acoustic far-field regimeallows for a much higher experimental flexibility,although at the cost of the straight forward in-terpretation of the measurements. More precisely,in the far-field, when the distance between detec-tor and source is large compared to the lateral ex-tend of the source, OA signals exhibit a diffraction-transformation which is characteristic for the un-derlying system parameters [14, 15, 16, 17]. In par-ticular, in the acoustic far-field, OA signals exhibita train of compression and rarefaction peaks andphases, signaling a sudden change of the absorptivecharacteristics of the underlying layered structure.In the presented article we thoroughly prepareand analyze polyvinyl alcohol hydrogel (PVA-H)phantoms, comprised of layers doped with differ-ent concentrations of melanin. The acoustic prop-erties of the PVA-Hs match those of soft tissue,i.e. human skin [18, 19]. Note that melanin is themain endogenous absorber in the epidermis [20],and, more importantly, in melanomas. Layers withhigher concentrations of melanin absorb a greateramount photothermal energy and expand more in-tensely than surrounding layers with low concen-trations. The stress waves emitted by these OAsources, detected in the acoustic far field after expe-riencing a shape transformation due to diffraction,are put under scrutiny here. Therefore, experimen-tal measurements are complemented by custom nu-merical simulations. Besides analytic theory andexperiment, the latter form a “third pillar” of con-temporary optoacoustic studies [5].The paper is organized as follows. In Sec. 2 we re-cap the theoretical background of optoacoustic sig-nal generation and detail our numerical approach tocompute the respective signals for point-detectors.In Sec. 3 we describe our experimental setup andelaborate on the preparation of the tissue phantomsused for our measurements, followed by details ofthe experiments and complementing simulations inin Sec. 4. We summarize our findings in Sec. 5.
2. Theory and numerical implementation
We briefly recap the general theory of optoacous-tic signal generation in Subsec. 2.1. Subsequently,in Subsec. 2.2, we customize the general optoacous-tic poisson integral to properly represent the lay-ered tissue phantoms and irradiation source profileused in our experiments. Finally, in Subsec. 2.3,we emphasize some important implications of the problem-inherent symmetries on our numerical im-plementation.
In thermal confinement, i.e. considering shortlaser pulses with pulse duration significantlysmaller than the thermal relaxation time of the un-derlying material [8, 21], the inhomogeneous optoa-coustic wave equation relating the scalar pressurefield p ( (cid:126)r, t ) to a heat absorption field H ( (cid:126)r, t ) reads (cid:104) ∂ t − c (cid:126) ∇ (cid:105) p ( (cid:126)r, t ) = ∂ t Γ H ( (cid:126)r, t ) . (1)Therein, c signifies the speed of sound and Γ refersto the Gr¨uneisen parameter, an effective parame-ter summarizing various macroscopic material con-stants, describing the fraction of absorbed heat thatis actually converted to acoustic pressure. As ev-ident from Eq. (1), temporal changes of the lo-cal heat absorption field serve as sources for stresswaves that form the optoacoustic signal. Followingthe common framework of stress confinement [5], weconsider a product ansatz for the heating functionin the form H ( (cid:126)r, t ) = W ( (cid:126)r ) δ ( t ) , (2)where W ( (cid:126)r ) represents the volumetric energy den-sity deposited in the irradiated region due to pho-tothermal heating by a laser pulse [22], which, onthe scale of typical acoustic propagation times, isassumed short enough to be represented by a delta-function. Consequently, an analytic solution for theoptoacoustic pressure at the field point (cid:126)r can be ob-tained from the corresponding Greens-function infree space, yielding the optoacoustic Poisson inte-gral [23, 24, 25] p ( (cid:126)r, t ) = Γ4 πc ∂ t (cid:90) V W ( (cid:126)r (cid:48) ) | (cid:126)r − (cid:126)r (cid:48) | δ ( | (cid:126)r − (cid:126)r (cid:48) | − ct ) d (cid:126)r (cid:48) , (3)where V denotes the “source volume” beyond which W ( (cid:126)r (cid:48) ) = 0 [26], and δ ( · ) limiting the integration toa time-dependent surface constraint by | (cid:126)r − (cid:126)r (cid:48) | = ct . As pointed out earlier, we consider non-scatteringcompounds, composed of (possibly) multiple plane-parallel layers, stacked along the z -direction of anassociated coordinate system. Whereas the acous-tic properties are assumed to be constant within2he medium, the optical properties within the ab-sorbing medium might change from layer to layer.Thus, the volumetric energy density naturally fac-tors according to W ( (cid:126)r ) = f f ( x, y ) g ( z ) , (4)wherein f signifies the energy fluence of the inci-dent laser beam on the z = 0 surface of the ab-sorbing material, and f ( x, y ) and g ( z ) specify thetwo-dimensional (2D) irradiation source profile andthe 1D axial absorption depth profile, respectively.Bearing in mind that we consider non-scatteringmedia, the latter follows Beer-Lambert’s law, i.e. g ( z ) = µ a ( z ) exp (cid:110) − (cid:90) z µ a ( z (cid:48) ) d z (cid:48) (cid:111) , (5)wherein µ a ( z ) denotes the depth-dependent absorp-tion coefficient.Note that, for a plane-normal irradiation sourcewith an axial symmetry, there are two useful aux-iliary reference frames based on cylindrical polarcoordinates: (i) Σ I where (cid:126)r = (cid:126)r ( ρ, φ, z ) with originon the beam axis at the surface of the absorbingmedium, and, (ii) Σ D where (cid:126)r (cid:48) = (cid:126)r (cid:48) ( ρ (cid:48) , φ (cid:48) , z (cid:48) ) withorigin at the detection point (cid:126)r D = ( x D , , z D ) in Σ I ,see Fig. 1. Both reference frames are related by thepoint transformation (cid:126)r (cid:48) ( ρ (cid:48) , φ (cid:48) , z (cid:48) ) = (cid:126)r − (cid:126)r D [27].In Σ I the irradiation source profile takes theconvenient form f ( x ( ρ, φ ) , y ( ρ, φ )) ≡ f I ( ρ ), wherebeam-profiling measurements for our experimentalsetup are consistent with a top-hat shape f I ( ρ ) = (cid:40) , if ρ ≤ a exp {− ( ρ − a ) /d } , if ρ > a . (6)In Σ D the constituents of the volumet-ric energy density read f D ( ρ (cid:48) , φ (cid:48) ) ≡ f ( x D + ρ (cid:48) cos( φ (cid:48) ) , ρ (cid:48) sin( φ (cid:48) )) and g D ( z ) ≡ g ( z (cid:48) − z D ), so thatthe optoacoustic Poisson integral, i.e. Eq. (3), takesthe form p D ( t ) = f Γ4 πc ∂ t (cid:90) (cid:90) (cid:90) V ρ f D ( ρ, φ ) g D ( z )( ρ + z ) / × δ (( ρ + z ) / − ct ) d ρ d φ d z. (7)Albeit the non-canonical formulation of the Poissonintegral in cylindrical polar coordinates might seema bit counterintuitive at first, it paves the way foran efficient numerical algorithm for the calculationof optoacoustic signals for layered media. Figure 1: Illustration of the two reference frames Σ I , withorigin on the beam axis at the surface of the absorbingmedium, and, Σ D with origin right at the detection point.Both coordinate systems are related by the transformation (cid:126)r (cid:48) ( ρ (cid:48) , φ (cid:48) , z (cid:48) ) = (cid:126)r − (cid:126)r D [27]. Considering cylindrical polarcoordinates in Σ D allows to factor the volumetric energydensity W ( (cid:126)r ) within the source volume V as detailed in thetext and to pre-compute the contribution of the irradiationsource profile along closed polar curves C ( ρ (cid:48) ) with radius ρ (cid:48) .This in turn yields an efficient numerical scheme to computethe optoacoustic signal p D ( t ) at the detection point (cid:126)r D . Considering a partitioningof the radial coordinate into N ρ equal sized values∆ ρ = L ρ /N ρ so that ρ i = i ∆ ρ with i = 0 . . . N ρ − W ( (cid:126)r ) in Σ D allows to pre-compute the con-tribution of the irradiation source profile in Eq. (7)along closed polar curves C ( ρ i ) with radius ρ i ac-cording to F D ( ρ i ) = lim N φ →∞ ρ i N φ − (cid:88) j =0 f D ( ρ i , φ j ) ∆ φ, (8)where ∆ φ = 2 π/N φ and φ j = j ∆ φ with j =0 . . . N φ −
1, thus completing the integration overthe azimuthal angle and providing the results in atabulated manner with time complexity O ( N ρ N φ ).This in turn yields an efficient numerical schemeto compute the optoacoustic signal p D ( t ) at thedetection point (cid:126)r D since the pending integrationscan, in a discretized setting with ∆ z = L z /N z sothat z k = k ∆ z for k = 0 . . . N z −
1, be carriedout with time complexity O ( N ρ N z ). Consequently,interpreting the δ -distribution in Eq. (7) as an in-dicator function that bins the values of the inte-grand according to the propagation time of the as-sociated stress waves, the overall algorithm com-pletes in time O ( N ρ N φ + N ρ N z ). Note that for3he special case x D = y D = 0, i.e. for detectionpoints on the beam axis, Eq. (8) further simplifiesto F D ( ρ i ) = 2 πρ i f I ( ρ i ), reducing the time complex-ity to only O ( N ρ N z ) [28]. During our numericalsimulations , for practical purposes and since weare only interested in the general shape of the op-toacoustic signal in order to compare them to thetransducer response, we set the value of the con-stants in Eq. (7) to f Γ /c ≡ π . Thus, the resultingsignal is obtained in arbitrary units, subsequentlyabbreviated as [a . u . ], making it necessary to adjustthe amplitude of the signal if we intent to comparethe results to actual measurements. Further, tomimic the finite thickness ∆ w of the transducer foil,see Sec. 3, we averaged the optoacoustic signal atthe detection point over a time interval ∆ t = ∆ w/c . Exemplary optoacoustic signals.
So as to facili-tate intuition and to display the equivalence ofthe numerical schemes implemented according toEqs. (3) and (7) in both, the acoustic near field(NF) and far field (FF), we illustrate typicaloptoacoustic signals in Fig. 2. Therein, the“cartesian” solver was based on a voxelized cu-bic representation of the source volume with side-lengths ( L x , L y , L z ) = (0 . , . , .
15) [cm] us-ing ( N x , N y , N z ) = (1500 , , L ρ , L z ) = (0 . , .
15) cm and ( N ρ , N φ , N z ) =(6000 , , a = 0 .
15 cm and d = a/
4. As finite thickness of the transducer foilwe considered ∆ w = 50 µ m in both setups.The dimensionless diffraction parameter D =2 | z D | / ( µa ) [14, 12] can be used to distinguish theacoustic near field (NF) at D <
D >
1. Here, we consider the effective parame-ters µ = (cid:104) µ a ( z ) (cid:105) and a = 1 . · a in case of multi-layered tissue phantoms. The simulations were per-formed at detection points on the beam axis, realiz-ing NF conditions with D ≈ .
15 at z D = − .
04 cmand FF conditions with D ≈ . z D = − . cτ = 0 . − . A Python implementation of our code for the solutionof the photoacoustic Poisson equation in cylindrical polarcoordinates, i.e. Eq. (7), can be found at [29]. c τ [cm] p ( t ) [ a . u . ] p ( t ) [ a . u . ] (a)(b) p D ( t ) - NF p D ( t ) - FF p (x D ,y D ,z D , τ ) 0 0.05 0.1 0.15 Figure 2: (Color online) Comparison of different solversfor the optoacoustic problem for layered media. The datacurves labeled by p D ( t ) refer to an implementation in cylin-drical polar coordinates according to Eq. 7. The curves arecomputed for a field point in the acoustic near-field (NF;red line) and far-field (FF; blue line) at z D = − .
04 cm and z D = − . p ( (cid:126)r D , τ )(black dashed lines). (a) Setup where the source-volume en-closes two absorbing layers consisting of µ a = 10 cm − inthe range z = 0 . − .
05 cm (light-gray shaded region) fol-lowed by µ a = 20 cm − in the range z = 0 . − . wave and accurately tracing the profile of the volu-metric energy density along the beam axis, followedby a pronounced diffraction valley for cτ > .
11 cm.The particular shape of the latter is characteristicfor the top-hat irradiation source profile used for thenumeric experiments. In contrast to this, as can beseen from Fig. 2, the FF signal features a successionof compression and rarefaction phases. Therein asudden increase (decrease) of the absorption coeffi-cient is signaled by a compression peak (rarefactiondip), cf. the sequence of peaks and dips at the points cτ = 0 , . , .
10 cm in Figs. 2(a,b). Further,the diffraction valley has caught up, forming rathershallow rarefaction phases in between the peaks anddips [30]. Finally, note the excellent agreement ofthe signals obtained by the two independent OAforward solvers.In a second series of simulations we clarified theinfluence of the radial deviation of the detectionpoint (cid:126)r D from the beam axis. Therefore we com-puted the excess pressure p D ( t ) at different posi-tions x D (cid:54) = 0 perceived in Σ I . The results for x D = 0 . x D = 0 . /e -width of the beam intensity profile, are il-lustrated in Fig. 3. As evident from Fig. 3(a), for4 D = − . D = 0 .
76 inthe acoustic NF, the optoacoustic signal appears tobe quite sensitive to the precise choice of x D . I.e.,as soon as the border of the plane-wave part of thesignal is approached, the transformation of the sig-nal due to diffraction is strongly visible. Comparingthe points z D = − . D = 11 .
4) in the “early”FF and z D = − . D = 19 .
0) in the “deep”FF, it is apparent that the optoacoustic signal inthe FF is less influenced by the off-axis deviationof the detection point, see Figs. 3(b,c). Also, notethat with increasing distance | z D | , the interjacentrarefaction phases level off and move closer to theleading compression peaks [30]. From the above weconclude that, if we complement actual measure-ments recorded in the FF via numerical simulations,we should find a good agreement between detectedand calculated signals even though both exhibit dif-ferent degrees of deviation from the beam axis.This completes the discussion of optoacoustic sig-nals and their generation from a point of view ofcomputational theoretical physics. Details regard-ing the optoacoustic detection device and the tissuephantoms are given in the subsequent section.
3. Methods and Material
Photoacoustic measurement setup.
In the follow-ing, the experimental setup is presented with focuson the phantom preparation process and arrange-ment of the layered tissue samples [31]. For thedetection of the OA pressure transient a self-builtpiezoelectric transducer is employed. This ultra-sound detector comprises a 9 µ m thick piezoelectricpolyvinylidenfluorid (PVDF) film on both sides ofwhich ∼
50 nm indium tin oxide (ITO) electrodesare sputtered [11]. The active area of the detectoris circular with a diameter of 1 mm. As acousticbacking layer a piece of hydrogel was prepared andplaced on top of the detector with a drop of distilledwater to ensure acoustic coupling. Due to identi-cal acoustic impedances of the backing layer andthe phantom in addition to the marginal extent ofthe PVDF film in comparison to the acoustic wave-lengths, the detector can be seen as acousticallytransparent [10].In contrast to the numerical approach followedin the preceding section, the irradiation in our ex-perimental setup had to be adjusted to an angleapproximately 20 ◦ off the plane normal, with thelight entering the phantom in close proximity ofthe detector, see Fig.4. As laser source, an optical c τ [cm] p D ( t ) [ a . u . ] p D ( t ) [ a . u . ] p D ( t ) [ a . u . ] (a)(b)(c) x D = 0.00.10.2 z D = -5.0 D = 19.0 z D = -1.0 D = 11.4 z D = -0.2 D = 0.76 x D = 0.00.10.2 z D = -5.0 D = 19.0 z D = -1.0 D = 11.4 0 0.05 0.1 x D = 0.00.10.2 z D = -5.0 D = 19.0 Figure 3: (Color online) Sensitivity of the optoacoustic sig-nal p D ( t ) on a radial deviation of the detection point (cid:126)r D fromthe beam axis, realized by setting x D (cid:54) = 0 cm, as explainedin the text. The subfigures refer to different distances z D ,where (a) z D = − . D = 0 .
76, (b) z D = − . D = 11 .
4) in the “early” FF,and, (c) z D = − . D = 19 .
0) in the “deep” FF. parametric oscillator (NL303G + PG122UV, Ek-spla, Lithuania) at a wavelength of 532 nm is cou-pled into a 800 µ m fiber (Ceramoptec, Optran WF800/880N). The pulse duration from the pump is3-6 ns. The beam profile measured after the fiberis in good agreement with a top-hat shape, whichis in accordance with the irradiation source profileEq. (6), considered for the previous numerical ex-periments, and parameters as detailed for the nu-merical simulations in Sec. 4.To improve the signal-to-noise ratio and matchthe electrical impedance, a custom build electri-cal pre-amplifier is connected to the detector elec-trodes. The voltages, corresponding to the detectedpressure, are recorded at 2 GSPS (Giga sample persecond) by a high-speed data acquisition card (Ag-ilent U1065A, up to 8GSPS). At such samplingrates, the expected ultrasound pressure profile ishighly over-sampled, thus, the point to point noisecan be smoothed out without loss of information. Aconservative estimation of the fastest changing sig-nal features yield a time window of 20 ns over whichsmoothing might be carried out, corresponding to40 consecutive data points. Polyvinyl alcohol based hydrogel tissue phantomrecipe.
The tissue phantoms used in our studies arecompounds composed of stacked layers of polyvinyl5lcohol hydrogel (PVA-H)[32]. The incentive to uti-lize PVA-H is its acoustic similarity to soft tissue,i.e. human skin [18]. In contrast to liquid phantomssuch as water ink solutions [12], hydrogels have theadvantage of being stackable without the need ofcontaining walls. Furthermore, while liquids wouldintermix at interfaces and thus require solid bound-aries in between, hydrogels allow sharp junctionsonly softened by diffusion. In the remainder thephantom creation recipe is detailed.Here, PVA-Hs are produced by mixing polyvinylalcohol granulate (Sigma-Aldrich 363146, Mw 85-124 99+% hydrolyzed) with distilled water at amixing ratio 1:5. Using a magnetic stirrer withheating, the dispersion is kept at 94 ◦ C for atleast 40 minutes while the stirring bar rotates at350 RPM, until it becomes a homogeneous solution.This viscous mass can be poured into any mold toobtain the desired form. Depending on the requiredthickness, a commercial metal spacing washer or a3D printed plastic ring of specific height was placedin between two glass plates, thus creating very flatPVA-H cylinders. Due to the much larger lateralextend of the phantoms, compared to the depth ofthe absorbing layers, boundary effects do not inter-fere with the optoacoustic signal.To facilitate polymerization the phantom is sub-jected to one freezing and thawing cycle . The phan-tom is placed in the freezer at − ◦ C for 2 days.Thawing is achieved by keeping the samples atroom temperature for a few minutes, afterwardsthe phantoms are ready to use. Crystallites pro-duced by freezing of water in the hydrogels wouldyield turbidity [33]. As a remedy, so as to obtainclear PVA-H, water soluble anti freezing agents areadded. Here, when the PVA is completely dis-solved, ∼
45 vol% pure ethanol was added to theaqueous solution incrementally, each time waitingfor the schlieren to dissolve. Especially after addingthe ethanol it is very important to keep the vesselclosed whenever possible.The optical properties of the samples can bemanipulated by inclusion of scatterers and or ab-sorbers. In our studies, synthetic melanin (Sigma-Aldrich, M0418-100MG) was chosen as absorber tomimic melanoma, that is, black skin cancer. Dueto its robustness to temperature, the finely groundmelanin can be included in the beginning of thephantom creation process, at the same time withthe PVA pellets. As a rough estimate it is as-sumed that melanomas contain as much melanin asAfrican skin. According to [34], dark skin contains ...
PI PII PIII (a)(b)
PVA hydrogelBacking layer ...
Figure 4: (Color online) Sketch of the experimental setup.(a) Arrangement of the components as discussed in the text,see sec. 3. (b) Layer composition of the three different phan-toms PI, PII and PIII used for our measurements and thenumerical simulations reported in sec. 4. The label “C” rep-resents clear PVA-H, “S” labels low absorption, and “M”stands for high absorption. Note that, the clear layer atthe bottom is 10 mm thick. Thus, signal reflections from in-terfaces of materials with differing acoustic properties occurwell outside the measurement range in that direction.
10 times the melanin as compared to very fair skin,not taking into account the type of melanin. So asto reproduce the contrast of a melanoma in Cau-casian skin the following different types of PVA-Hlayers were created:(i) PVA-H without melanin, referred to as “C”,(ii) PVA-H with 1 mg/mL of melanin, referred toas “M”, and,(iii) PVA-H with 0 . Further hints for handling the tissue phantoms.
Be-low we list practical hints and findings which proved6ery helpful in the course of tissue phantom produc-tion:H1: In a closed screw neck bottle the aqueous solu-tion can be stored for weeks at room temper-ature. However, note that after the flask hasbeen reheated and opened several times, in-evitable ethanol evaporation might cause tur-bidity upon freezing.H2: The undesirable formation of bubbles withinthe phantoms can be avoided, to a greatextend, by overfilling the spacing ring andmounting the top glass plate on the sample af-ter the hydrogel has settled for a while. Thisprocedure prohibits trapping of air at the spac-ing ring boundaries as well as giving the hydro-gel some time to degas.H3: While stacking the PVA-H layers in prepara-tion for a measurement, the individual phan-tom layers should be kept wet by means of dis-tilled water in order to prevent them from sick-ing together with one another and, most of all,themselves. Also, a proper watery film pro-hibits the inclusion of air in between layers.
4. Results
Below we complement measured PA signals, ob-tained from measurements on the three tissue phan-toms PI − III, discussed in Sec. 3 and illustratedin Fig. 4(b), with custom simulations obtained interms of the numerical framework detailed in Sec.2. As evident from the comparison of the experi-mental setup with the simulation framework, thereare three distinctions between experiment and the-ory which have to be kept in mind while interpret-ing the results: (i) while the irradiation source isassumed to be plain normal incident for our sim-ulations, the direction of incidence in the experi-mental setup exhibits a nonzero angle off the planenormal. Additionally, due to unavoidable refractionat the phantom surface, the beam profile is likelyto be non-symmetric and slightly divergent. Hence,the top-hat beam shape assumed in our simulationapproach can only been seen as an approximationof the experimental conditions. (ii) although it isprobable that all the measurements are performed,at least to some extend, off-axis we opt for modelingand numerical simulations in an on-axis approach.As demonstrated in Subsec. 2.3 and illustrated inFig. 3(c), we expect the principal signal shape in c τ [cm] z [cm] X ( t ) [ a . u . ] X ( t ) [ a . u . ] X ( t ) [ a . u . ] (a)(b)(c) X ( t ) = V T ( t ) p D ( t )(a)(b)(c)-0.1 0 0.1 0.2 0.3(a)(b)(c) Figure 5: (Color online) Comparison of optoacoustic sig-nals for layered media obtained from measurements ( labeled“ V D ( t )”; orange solid lines) and numerical simulations (la-beled “ p D ( t )”; blue dashed lines). The top and bottom ab-scissas refer to the distance z traveled by the signal and theretarded signal depth cτ = ct − z D (where z D = 0 . the acoustic far-field to change only at a small rateupon deviation from the beam axis. (iii) as pointedout in sec. 3, the active area of the transducer has aradius of 0 . Comparison of optoacoustic signals obtained in the-ory and experiment.
The measured optoacousticsignals for the tissue phantoms PI − III along withthe simulated curves are illustrated in Figs. 5(a-c). In principle all three measured curves exhibitthe characteristic features expected for signals ob-served in the acoustic far-field. Thus these mea-surements are well suited for the purpose of optoa-coustic depth profiling [12]. In particular, for the7imulation of PI we considered a single layer withabsorption coefficient µ a = 11 cm − in the range z = 0 . − .
395 cm (note that z is measured withrespect to the origin of Σ D ), indicated by a grayshaded region representing a highly absorbing layer(introduced as “M” in sect. 3). The top-hat beamshape parameters within the simulation where setto a = 0 .
054 cm and R ≡ d/a = 1 .
5. Note that bothprincipal features of the signal, i.e. the initial com-pression peak as well as the trailing rarefaction dipare reproduced well and the intermediate rarefac-tion phase matches well in theory and experiment.The subsequent long and shallow rarefaction phasefor z > . µ a = 1 . − , i.e. type-S, in the range z = 0 . − .
408 cm (light-gray shaded region), fol-lowed by a type-M layer with µ a = 11 cm − in therange z = 0 . − .
504 cm (gray shaded region).Therein, the beam shape parameters where set to a = 0 .
056 cm and R = 1 .
2. Here, all three ex-pected characteristic signal features, i.e. the initialsmall compression peak, the interjacent high com-pression peak as well as the trailing rarefaction dipmatch well for theory and experiment.Finally, phantom PIII was modeled by consid-ering a type-S layer with µ a = 1 . − in therange z = 0 . − . µ a = 11 cm − inthe range z = 0 . − .
595 cm (gray shaded region).Therein, the beam shape parameters where set to a = 0 .
08 cm and R = 1 .
2. Again, all three charac-teristic signal features are reproduced well by the-ory and experiment.Note that, as pointed out in subsec. 2.3, it is nec-essary to adjust the scale of the amplitude of thecomputed PA signal if we intend to compare it tothe transducer response. The respective scaling fac-tor was obtained from the simulated and measuredcurves for tissue phantom PI and was subsequentlyused in the other two cases to achieve the excellentagreement displayed in Figs. 5(a-c).
Reconstruction of the initial volumetric energy dis-tribution.
As discussed in the literature, in theacoustic far-field, the observed optoacoustic signalcan be related to the initial volumetric energy den-
0 0.1 0.2 0.3 0.4 0.3 0.4 0.5 0.6 0.7 c τ [cm] z [cm] p , FF ( c τ ) [ a . u . ] (a)(b) T - p (c τ ) E - z D = -0.3 cmT - z D = -0.3 cmT - z D = -0.9 cmT - z D = -4.0 cm10 -8 -7 -6 z D [cm]MSE(a)(b) Figure 6: (Color online) Reconstruction of the initial volu-metric energy density in the far-field (FF) approximation fortissue phantom PIII, see Fig. 5(c). (a) comparison of the ex-act initial distribution of acoustic stress p (black solid line;labeled “T”) to FF reconstructed predictors p , FF simulatedat different measurement points z D (blue lines; labeled “T”)and the FF reconstructed predictor derived from our mea-surement (orange solid line; labeled “E”). (b) mean squareerror MSE between the exact initial volumetric energy den-sity and the FF reconstructed value from the optoacousticsignals calculated at different detection points at z D = − . − sity by means of a temporal derivative [12, 14].Consequently, as discussed in Ref. [12], this of-fers the possibility to reconstruct the initial acous-tic stress distribution p ( z ) = Γ W ( z ) in the limit D (cid:29)
1. Albeit Ref. [12, 13] used the integral ofthe measured acoustic signals as a visual aid forimaging purposes, cf. Fig. 9(c) of Ref. [12], and Fig.8(b) of Ref. [13], they did not elaborate on this is-sue any further. Albeit we agree that the FF sig-nals are naturally suited for the purpose of optoa-coustic depth profiling, we here attempt to explorethe use of the above idea in order to obtain a pre-dictor p , FF ≈ p in terms of a FF approximationfor tissue phantom PIII. This is illustrated in Fig.6(a), where we show the exact initial distributionof acoustic stress p (solid black line) by means ofwhich the numerical simulations where carried out,together with the FF reconstructed predictors p , FF simulated at three different measurement points z D = − . , − . , − . z D = − . | z D | yieldsa more consistent estimate. In the limit | z D | → ∞ this is limited only by the temporal averaging of thesignal, implemented to mimic a finite thickness ofthe transducer foil.This can be assessed on a quantitative basisby monitoring the mean squared error MSE = (cid:80) N z − i =0 [ p ( z i ) − p , FF ( z i )] /N z in a discretized set-ting, with z i as in Subsec. 2.3, see Fig. 6(b). Notethat in advance, the above signals are normalized inorder to ensure (cid:80) i X ( z i ) = 1 for both, X = p and p , FF . As evident from the figure, the MSE mightbe reduced by a solid order of magnitude upon mov-ing the signal detection from z D = − . − .
5. Summary and Conclusions
In the presented article we discussed an efficientnumerical procedure for the calculation of optoa-coustic signals in layered media, based on a nu-merical integration of the optoacoustic Poisson inte-gral in cylindrical polar coordinates, in combinationwith experimental measurements on PVA basedhydrogel tissue phantoms. In summary, we ob-served that far-field measurements on tissue phan-toms composed of layers with different concentra-tions of melanin are in striking agreement with cus-tom numerical simulations and exhibit all the char-acteristic features that allow for optoacoustic depthprofiling. Further, in our experiments, the signal tonoise ratio of s ingle measurements was sufficientlyhigh to omit any signal post-processing. In con-trast to the experimental measurements, the simu-lations are performed with on axis illumination andassuming an ideal pointlike detector. Nonetheless,simulation and experiment agree very well over all,which highlights the robustness of the signal anal-ysis and simulation against small deviations. Fi-nally, we showcased the possibility to reconstructthe initial pressure profile in a far-field approxima-tion by numerical integration. Even though exactreconstruction would require an ideal detector inaddition to an infinite distance between source anddetector, the pressure profile reconstructed here (atfinite distance | z D | = 1 cm and finite detector radius0 . Acknowledgments
We thank J. Stritzel for valuable discussionsand comments, as well as for critically readingthe manuscript. We further thank M. Wilke forassisting in the preparation of the PVA-H tissuephantoms. E. B. acknowledges support from theGerman Federal Ministry of Education and Re-search (BMBF) in the framework of the projectMeDiOO (Grant FKZ 03V0826). O. M. ac-knowledges support from the VolkswagenStiftungwithin the “Nieders¨achsisches Vorab” program inthe framework of the project “Hybrid NumericalOptics” (Grant ZN 3061). Further valuable discus-sions within the collaboration of projects MeDiOOand HYMNOS at HOT are gratefully acknowl-edged.
References [1] L. V. Wang, S. Hu, Photoacoustic tomography: in vivoimaging from organelles to organs, Science 335 (6075)(2012) 1458–1462.[2] I. Stoffels, S. Morscher, I. Helfrich, U. Hillen, J. Leyh,N. C. Burton, T. C. P. Sardella, J. Claussen, T. D.Poeppel, H. S. Bachmann, A. Roesch, K. Griewank,D. Schadendorf, M. Gunzer, J. Klode, Metastatic sta-tus of sentinel lymph nodes in melanoma determinednoninvasively with multispectral optoacoustic imaging,Science Translational Medicine 7 (317) (2015) 317ra199.[3] I. Stoffels, S. Morscher, I. Helfrich, U. Hillen, J. Leyh,N. C. Burton, T. C. P. Sardella, J. Claussen, T. D.Poeppel, H. S. Bachmann, A. Roesch, K. Griewank,D. Schadendorf, M. Gunzer, J. Klode, Erratum forthe research article: “metastatic status of sentinellymph nodes in melanoma determined noninvasivelywith multispectral optoacoustic imaging” by i. stoffels,s. morscher, i. helfrich, u. hillen, j. lehy, n. c. burton,t. c. p. sardella, j. c. . . , Science Translational Medicine7 (319) (2015) 319er8.[4] J. Gateau, A. Chekkoury, V. Ntziachristos, Ultra-wideband three-dimensional optoacoustic tomography,Optics Letters 38 (22) (2013) 4671.[5] L. Wang, Photoacoustic Imaging and Spectroscopy, Op-tical Science and Engineering, CRC Press, 2009.[6] D. Colton, R. Kress, Inverse Acoustic and Electromag-netic Scattering Theory (3rd Ed.), Springer, 2013.[7] P. Kuchment, L. Kunyansky, Mathematics of ther-moacoustic tomography, European Journal of AppliedMathematics 19 (2008) 191–224.
8] R. A. Kruger, P. Liu, Y. Fang, C. R. Appledorn, Pho-toacoustic ultrasound (paus) - reconstruction tomogra-phy, Medical Physics 22 (1995) 1605–1609.[9] H. F. Zhang, K. Maslov, G. Stoica, L. V. Wang, Func-tional photoacoustic microscopy for high-resolution andnoninvasive in vivo imaging, Nature Biotechnology24 (7) (2006) 848–851.[10] M. Jaeger, J. J. Niederhauser, M. Hejazi, M. Frenz,Diffraction-free acoustic detection for optoacousticdepth profiling of tissue using an optically transparentpolyvinylidene fluoride pressure transducer operated inbackward and forward mode, Journal of Biomedical Op-tics 10 (2) (2005) 024035.[11] J. J. Niederhauser, M. Jaeger, M. Hejazi, H. Keppner,M. Frenz, Transparent ito coated pvdf transducer foroptoacoustic depth profiling, Optics Communications253 (4-6) (2005) 401–406.[12] G. Paltauf, H. Schmidt-Kloiber, Pulsed optoacousticcharacterization of layered media, Journal of AppliedPhysics 88 (2000) 1624–1631.[13] G. Paltauf, H. Schmidt-Kloiber, Optoakustische spek-troskopie und bildgebung, Zeitschrift f¨ur MedizinischePhysik 12 (2002) 35–42.[14] A. Karabutov, N. B. Podymova, V. S. Letokhov, Time-resolved laser optoacoustic tomography of inhomoge-neous media, Appl. Phys. B 63 (1998) 545–563.[15] M. W. Sigrist, Laser generation of acoustic waves inliquids and gases, J. Appl. Phys. 60 (1986) R83–R122.[16] G. J. Diebold, M. I. Khan, S. M. Park, PhotoacousticSignatures of Particulate Matter: Optical Productionof Acoustic Monopole Radiation, Science 250 (4977)(1990) 101–104.[17] G. J. Diebold, T. Sun, M. I. Khan, Photoacousticmonopole radiation in one, two, and three dimensions,Phys. Rev. Lett. 67 (1991) 3384–3387.[18] A. Kharine, S. Manohar, R. Seeton, R. G. M. Kolkman,R. A. Bolt, W. Steenbergen, F. F. M. deMul, Poly(vinylalcohol) gels for use as tissue phantoms in photoacousticmammography, Physics in Medicine and Biology 48 (3)(2003) 357–370.[19] C. M. Moran, N. L. Bush, J. C. Bamber, Ultrasonicpropagation properties of excised human skin, Ultra-sound in Medicine & Biology 21 (9) (1995) 1177–1190.[20] J. E. Smit, A. F. Grobler, R. W. Sparrow, Influenceof variation in eumelanin content on absorbance spec-tra of liquid skin-like phantoms, Photochemistry andphotobiology 87 (1) (2011) 64–71.[21] Note that Ref. [8] put the heat conduction equation un-der scrutiny, finding that for pulse durations t p < µ sand absorption lengths (cid:96) > ≈ . Here, we consider absorptionlengths of (cid:96) ≈ − t p ≈