**1. Introduction**

One of the goals of modern medical imaging is to simultaneously provide molecular, functional, and structural/anatomical information corresponding to various tissues. To achieve such comprehensive information, a combination of conventional imaging technologies such as X-ray computed tomography (CT) [1], positron emission tomography (PET) [2] and magnetic resonance imaging (MRI) [3] are used. CT provides anatomical contrast, PET provides molecular and metabolic information, and MRI provides both functional and anatomical contrasts. A PET–CT or PET–MRI combination is therefore widely used for simultaneously mapping molecular and anatomical contrasts [4–6].

All of these methods have significant downsides: one of which is the financial burden, which means that these methods are not ideal for routine imaging, and a second is the ionizing radiation used in CT and PET. As an alternative, ultrasound imaging that uses non-ionizing radiation is routinely used in several clinical applications for anatomical imaging [7]. However, it lacks the molecular or functional information necessary for detecting the early symptoms of disease.

Alternative medical imaging modalities that provide anatomical, functional, and molecular information about the tissue, while being lower cost and not being as restrictive to patient movement, are needed. Photoacoustic computed tomography (PACT) and radio frequency (RF)-induced

acoustic computed tomography (RACT), together known under the general term thermoacoustic computed tomography (TACT), match these objectives of lower cost functional/molecular imaging, using non-ionizing electromagnetic radiation [8–16]. While PACT maps optical absorption contrast using optical radiation induced acoustic wave detection, RACT maps tissue conductivity using RF induced acoustic wave detection. More importantly, since these hybrid (combining electromagnetic radiation and acoustic detection) imaging modalities share the same ultrasound detection platform, combined trimodality PACT–RACT–UCT (ultrasound computed tomography) systems have been realized for mapping functional, molecular, and anatomical contrasts [17]. The molecular absorption of electromagnetic energy causes thermal expansion in the tissue, which then leads to generation of acoustic waves. The acoustic waves propagate out of the tissue and are received by ultrasound transducers located at the boundary of the body. This data is then used to reconstruct thermoacoustic images displaying electromagnetic absorption contrast at ultrasonic spatial resolution. The imaging depth and spatial resolution in TACT is scalable with the frequency of excitation radiation and ultrasound transducer. The fact that the detected acoustic signal arises directly from specific molecules inside the tissue makes TACT a molecular/functional imaging technology. In PACT, the tissue chromophores—such as oxy-hemoglobin, deoxy-hemoglobin, melanin, and lipids—absorb light photons in the wavelength range from 400 nm to 1200 nm. By using different wavelengths, which takes advantage of the resonance peaks in the absorption spectrum of the imaged molecules, the distribution of different molecules inside the tissue can be mapped. In RACT, a radio frequency source in the frequency range from 434 MHz to 9 GHz is used for mapping the water distribution. The thermoacoustic effect in this frequency range is dependent on the conductivity distribution of the medium. The conductivity difference between water, tissue, and tumors can then give useful functional images.

Accurate numerical modeling of TACT is paramount for the development of robust reconstruction algorithms to quantify the electromagnetic absorption properties of the tissue. In PACT, the forward optical simulation of total light fluence ( Φ), calculated inside the tissue medium, is usually achieved using either Monte Carlo simulations or the software package NIRFast, which solves the light diffusion equation. The fluence distribution is then converted to an initial pressure rise, which is further propagated through the tissue medium and detected by the ultrasound transducers located on the boundary using acoustic simulation tools such as the K-wave toolbox. NIRFast [18] uses the finite element method (FEM), while K-wave uses a spectral-based finite difference method (FDM) to model the propagation of acoustic waves [19]. The goal of the reconstruction problem is then to recover the tissue properties (optical absorption in PACT and conductivity in RACT) given only the sensor data.

Groups that use a single simulation platform for both the optical and acoustic propagation needed for the hybrid PACT/RACT technique are less common. Recently, there are studies that use the commercial software COMSOL to solve both propagation problems with a single software package [20,21]. As an alternative to the commercial software, the simulation system described in this paper uses the open source softwares Gmsh [22] and GetDP [23], often combined under the name ONELAB [24]. ONELAB is an FEM solver, which uses Gmsh for creating the FEM mesh, and GetDP for solving generic partial differential equations (PDEs) with the FEM method. Advantages of using Gmsh include its ability to create user defined meshes, but also having standard interfaces with other commonly used mesh and computer-aided design (CAD) software such as STEP, IGES, and STL. Segmented DICOM images commonly used in MRI can then be converted to a mesh that Gmsh understands, creating realistic phantoms on which to test algorithms. The ability of GetDP to solve generic PDEs allows the user to implement algorithms to solve for optical and acoustic propagation, similar to the combination of NIRFast and K-wave. Both propagation methods, as well as any reconstruction methods, are implemented in the same mesh, with no loss of precision. Since there is less chance of numerical error, a single software platform represents a more accurate approximation of a real-world scenario. Besides PACT with optical sources, we demonstrate that the generality of GetDP allows for sources in the RF regime of the electromagnetic spectrum to also be simulated

for RACT. The use of the ONELAB software package allows accurate TACT modeling, in order to develop algorithms for functional imaging of the human body. Although there are experimental studies combining PACT and RACT techniques into one setup, perhaps our study is the first to report a single TACT simulation platform for simulating both PACT and RACT.

The rest of the paper is organized as follows: In Section 2, we describe methods and materials. Sections 3.1 and 3.2 show PACT and RACT results on a homogeneous phantom. Sections 3.3 and 3.4 show PACT and RACT results on an approximate breast phantom. Section 4 is a discussion of the results in Section 3, comparing errors in reconstruction between PACT and RACT, and between the two phantoms. Section 5 concludes the paper.

#### **2. Materials and Methods**

There are three main phenomena that need to be modeled in thermoacoustics. The first is the propagation of the initial energy source, in this study near infrared optical or radio frequency electromagnetic waves, and calculation of the total fluence/intensity distribution inside the tissue medium. The second phenomenon is the acoustic propagation, after calculating an initial pressure from the intensity distribution. The third phenomenon is the reconstruction of the tissue parameters, which depend on the type of energy source used for tissue excitation. A flowchart showing the similar steps between RACT and PACT is shown in Figure 1. This study is focused on using a time reversal [25,26] reconstruction method to recover the initial pressure, followed by a simple division to find the tissue parameters. There are several alternative approaches for reconstructing the tissue parameters that can be integrated with our study in the future. Back-projection methods may be used when using specific geometries [27] and parametrix methods for more general geometries [28]. The back-projection methods are analytically exact, but only for geometries with a special symmetry, and with a constant speed of sound. Parametrix methods do not provide analytically exact methods, but instead give approximations with error bounds for general geometries, and can even be developed for regions with varying speeds of sound. A survey of several reconstruction methods is given in Reference [29].

**Figure 1.** The similarities between key steps of the photoacoustic computed tomography (PACT, blue) and radio frequency (RF)-induced acoustic computed tomography (RACT, red) workflows, with methods that are shared represented in yellow. Step A represents the finite element mesh in Section 2.5, used for both types of simulation. Step B1 is represented in Figure 9a, while step B2 is represented in Figure 7a. Step C2 uses the optical diffusion Equation (1). Step C1 uses the electromagnetic wave Equation (5). The resulting pressure rise for the optical (D2) and RF (D1) radiations are represented by Figures 7b and 9b, respectively. The initial pressure propagates out to the boundary (E) via the scalar wave Equation (7). The data received by acoustic sensors at the boundary (F) can then be used in reconstruction algorithms (G) to recover the initial pressure. In this paper, a time-reversal algorithm is applied for reconstruction (Section 2.4). The reconstructed pressure is then divided by the intensity (I2, I1) to recover the recover the absorption or conductivity for the optical and RF cases, shown in Figures 7c and 9c, respectively.

#### *2.1. Optical Propagation for Photoacoustic Computed Tomography (PACT)*

A full description of optical propagation would use the radiative transport equation in its full generality. This equation can be modeled by the use of Monte Carlo methods, with significant computational cost [30,31]. Since the the photoacoustic effect depends on the optical fluence rate, and the wavelengths used correspond to the near infrared (NIR) regime, a diffusion approximation to the radiative transport equation is often used [32]. The diffusion equation for the fluence rate is:

$$-\nabla \cdot \mathbf{x}(\mathbf{x}) \nabla \Phi(\mathbf{x}, \lambda) + \mu\_a(\mathbf{x}) \Phi(\mathbf{x}) = q(\mathbf{x}) \ \mathbf{x} \in \Omega \tag{1}$$

$$
\Phi(\mathbf{x}) = \mathbf{0} \; \mathbf{x} \in \mathfrak{Y} \Omega. \tag{2}
$$

In this equation, Φ is the optical fluence rate, or intensity of the light. *μa* is the absorption coefficient at each point in space, while *κ* = 1 <sup>3</sup>(*μa*+*μs*) is the diffusion coefficient, calculated using the absorption coefficient and the reduced scattering coefficient *μs*. *q* is the source function, the initial laser pulse that irradiates the tissue medium. Ω is the tissue domain of interest, usually a subset of R*<sup>n</sup>*, with *n* corresponding to the dimension of the problem (2 or 3). Equation (2) corresponds to the Dirichlet boundary condition used in this study. The amount of pressure generated by the optical energy is given by a constant of proportionality, the Gruneisen parameter Γ = *<sup>α</sup>v*2*s Cp* , where *α* is the volume thermal expansion coefficient, *vs* is the speed of sound, and *Cp* is the heat capacity at constant pressure [33]. Therefore, for a given fluence rate Φ, the initial pressure is given by:

$$
\mu\_0(\mathbf{x}) = \Gamma(\mathbf{x}) \Phi(\mathbf{x}) \mu\_a(\mathbf{x}). \tag{3}
$$

While Γ will physically vary over space, this small deviation is ignored in this study, and the parameter is assumed constant, with a value of 0.1, which approximates standard tissue.

#### *2.2. Radio Frequency Propagation for RF Acoustic Computed Tomography (RACT)*

Radio frequency radiation is modeled by Maxwell's equations, which when modeled at a single frequency can be reduced to the following wave equation:

$$
\nabla \times \nabla \times E(\mathbf{x}) - \omega^2 \mu(\mathbf{x}) \varepsilon(\mathbf{x}) E(\mathbf{x}) = q(\mathbf{x}) \cdot \mathbf{x} \in \Omega. \tag{4}
$$

In this equation, *E* is the complex electric field, *ω* = 2*π f* is the frequency, and *μ* and *ε* are the permeability and permittivity of the medium, respectively. In a material with electric loss, the permittivity is complex valued, with *ε* = *ε*real − *j σω* , and *σ* = 0 for a medium with no loss. When the simulation is restricted to two dimensions (assuming a *TMz* polarization), Equation (4) can be reduced to a scalar Helmholtz equation:

$$
\nabla^2 E\_3(\mathbf{x}) + k(\mathbf{x})^2 E\_3(\mathbf{x}) = q(\mathbf{x}) \cdot \mathbf{x} \in \Omega,\tag{5}
$$

where *k* = *ωc* is the wavenumber, and can be complex in a tissue with energy loss. For boundary conditions, a Dirichlet boundary condition similar to Equation (2) is implemented. However, this would not prevent the electromagnetic waves from reflecting off the boundaries of the domain, causing numerical artifacts. Therefore, while there is a Dirichlet boundary condition, there is also a non-physical space occupied by a perfectly matched layer (PML) to absorb the outgoing waves and prevent reflection [34]. Once the electric field is calculated, the initial pressure is given by three contributions [35]. The first contribution is due to the conductivity of the medium, *pcond* = *V σ*2 |*E*|<sup>2</sup>*dV*. The second due to the permittivity *V ε*0*εr* 2 |*E*|<sup>2</sup>*dV*, and the third the permeability *V μ*0*μr* 2 |*H*|<sup>2</sup>*dV*. In practice, for the tissue media of interest in thermoacoustic tomography, the contribution to the pressure due to the conductivity dominates over the permittivity and permeability, and so these terms are

generally ignored. Therefore, similar to the optical absorption in Equation (3), the initial pressure due to the electric field is given as:

$$p\_0(\mathbf{x}) = \Gamma(\mathbf{x}) \frac{\sigma(\mathbf{x})}{2} |E(\mathbf{x})|^2. \tag{6}$$
