Effect of HIFU-Induced Thermal Ablation in Numerical Breast Phantom

Abstract

Breast cancer is a leading cause of cancer-related deaths in women, and treatment involved invasive surgery such as lumpectomy. In the last decade, a non-invasive, non-contact high-intensity focused ultrasound (HIFU) therapy was developed for treatment with promising results. However, its success rate depends on patient selection, tissue heterogeneities, HIFU operational parameters, and even imaging techniques. In this emerging field, computer simulations can provide us with a much-needed platform to learn, test, and deduce results virtually before conducting experiments. In this study, we used three different classes of anatomically realistic numerical breast phantoms from clinical contrast-enhanced magnetic resonance imaging (MRI) data, including scattered-, heterogeneous-, and extremely dense-type breasts. Upon assigning the appropriate acoustic and optical parameters to the tissues within, we simulated HIFU propagation by using the k-Wave toolbox in MATLAB and compared the changes introduced in the three types of breasts. It was found that scattered-type breast was best-suited for HIFU therapy. Furthermore, we simulated light-beam propagation with the ValoMC toolbox in MATLAB after introducing the lesion to compare the distribution of the initial pressure generated via the photoacoustic effect. This simulation study will be of significant clinical impact, especially in the study and management of HIFU-based treatments, which are individual/tissue-selective in nature.

1. Introduction

Breast cancer is the leading cause of cancer in women. According to statistics, $12.5\%$ of all cancer diagnoses among women occur from this deadly disease [1]. Breast cancer occurs when essential genes mutate in the breast, leading to uncontrolled cell growth [2]. Several genetic and environmental features also influence such mutations and tumour mass growth. In this regard, early detection and efficient treatment can play a crucial role in reducing the mortality rate caused by this deadly disease. Presently, mammography and ultrasonography are routinely used in clinics as screening methods [3,4]. Both modalities deliver morphological information about the underlying tissue. Recently, breast-conserving surgery combined with adjuvant therapies has been the standard early-stage therapeutic approach, which has been shown to be equally effective as mastectomy [5,6]. Breast surgeons normally use oncoplastic technique to decrease the impact of the local tumour excision on cosmesis. Various factors have been encountered that produce negative cosmesis consequences, such as being overweight, having tiny or very large breasts, developing a tumour in the lower quadrant of the breast, and having a large breast tumour [7]. Many therapeutic modalities are utilised to treat cancerous breast tissue, such as chemotherapy, radiotherapy, and hormone therapy [8]. Sometimes, these methods increase the harmful surgical impact due to tissue injuries, such as shrinking, fibrosis, and breast deformity. These are also non-invasive in nature. A variety of minimally invasive ablation techniques, such as radiofrequency ablation, cryoablation, and microwave ablation, have been introduced as alternatives to the state-of-the-art techniques [9,10,11]. All these techniques have promising attributes with improved cosmetic outcomes and mortality. However, the percutaneous insertion of a probe is required for all these methods [12,13]. Therefore, there is a clinical need for an entirely non-invasive ablative treatment for the betterment of medical science.

High-intensity focused ultrasound (HIFU) is non-invasive thermal ablation that has been employed as an effective and safe treatment in many biomedical applications, such as bone metastases [14], uterine fibroma [15], liver tumours [16,17], and pancreatic tumours [18], among others. In HIFU therapy, a high-intensity ultrasound beam emitted from a specially designed (high-intensity) transducer is focused onto a particular tissue region [19,20]. The sharply focused acoustic energy is transformed into thermal energy, which leads to tissue coagulation or ablation in that specific tissue region without damaging the surrounding tissue [21]. These advantageous features are an excellent therapeutic pathway for breast cancer treatment. In past studies, it has been shown that, compared with lumpectomy, HIFU has potentially improved cosmetic outcomes, prevented scarring possibilities, and reduced postoperative recovery time. In addition, the risk of surgical complications is eliminated due to the non-invasive nature of the therapy. However, the shape, size, and location of the tissue ablation depend on the acoustic properties and structure of tissues within the breast, along with the acoustic parameters of the HIFU. The success rate of HIFU ablation also depends on tissue heterogeneities, imaging techniques, the ablation protocol, and patient selection [22].

HIFU-induced breast tissue coagulation is traditionally monitored by X-ray mammography [23]. However, it has decreased sensitivity in dense breast tissue. As an alternative, magnetic resonance imaging (MRI) has good sensitivity and specificity in dense breast tissue [24]. However, it carries the burden of high cost and processing time. Moreover, MRI requires an external contrast agent, which may lead to an allergic reaction in the central nervous system, as well as renal damage in some cases. Ultrasound imaging, on the other hand, is cost-effective and radiation-free but has low contrast and sensitivity towards the coagulated region [25]. Thus, a new imaging method that is not affected by breast density while overcoming the limitations mentioned above is demanded. Recently, photoacoustic (PA) imaging has attracted considerable attention in breast imaging applications [26,27]. It is an emerging modality that offers deep tissue imaging with rich optical contrast and scalable spatial resolution that is much higher than that of any purely optical imaging modalities [28,29]. Additionally, PA imaging enables the quantification of functional parameters (haemoglobin concentration and oxygen saturation) of breast tissue during cancer treatment [30].

Numerical simulation can be a useful tool to learn about all the above aspects. In the present study, we simulated acoustic and optical propagation inside a numerical breast phantom by employing the k-Wave [31] and ValoMC toolboxes [32] in MATLAB [33]. This simulation will be helpful for the optimal arrangement of HIFU with respect to an individual’s breast structural composition. Identification of suitable HIFU operational parameters via simulations performed here will greatly aid in the success of HIFU therapy. To this end, we used anatomically realistic numerical breast phantoms obtained from clinical contrast-enhanced magnetic resonance (MR) imaging data [34]. We then used the k-Wave toolbox to virtually simulate HIFU propagation. The simulations, being specific to a patient’s breast structure and composition, can be highly effective in deciding the operational parameters.

There have been many reported studies of HIFU simulation in breasts. Here, we extend this technique by comparing the effect of HIFU on three different types of breasts based on the density and distribution of fibroglandular tissue and discuss the potential consequences of the results obtained if carried out in the real world. We further simulated light-beam propagation after introducing lesions to compare initial pressure distributions due to the photoacoustic effect in the three breast types. This will be of great impact on photoacoustic imaging and the reconstruction of lesions in our future works.

2. Materials and Methods

2.1. Numerical Breast Phantom

The breast models (corresponding to free-hanging position) were reproduced from a publicly available online resource—the ‘Optical and Acoustic Breast Phantom Database (OA-Breast)’ [34,35]. Typically, these are contrast-enhanced MR data from healthy patients where contrast agents are injected to provide higher contrast for blood vessels and other structures. The MR data were then processed and segmented into four different tissue types: skin layer, blood vessels, fat, and fibroglandular tissue. Readers can refer to the sources cited above for details. Three sample types of breasts are classified as scattered fibroglandular, heterogeneously dense-, and extremely dense-type breasts. By assigning specific tissue optical and acoustic properties for each tissue type within the breast, we can generate a numerical breast phantom. The voxel size of all three sample breast phantoms is 0.2 mm. The loaded digital phantoms have a resolution of 0.2 mm as well. The 2D cross sections used for HIFU simulations have the same resolution in all three samples. In our study, a 2D model was used. Figure 1 shows a 2D cross-section of the extremely dense breast sample extracted from the 3D (mean) projected breast. The dimensions of the breast are 615 × 752 × 495 pixels (2D cross-section dimensions: 615 × 752 pixels, for an effective physical size of breasts of 12.3 cm × 15.04 cm). The cross-section was chosen roughly at the midpoint of the breast along the z-axis. Values used for acoustic, thermal, and optical properties in the simulation are shown in

Table 1. Acoustic parameters corresponding to various parts of breast tissues (skin, blood vessels, fat, and glands) and the surrounding water employed in simulation for HIFU wave propagation [36].

Characteristics	Skin	Vessel	Fat	Gland	Water
Attenuation coefficient [dB/(MHzy·cm)]	1.84	0.61	0.34	0.75	0.002
Power factor [y]	1.69	1.3	0.73	1.5	2
Sound speed [m/s]	1650	1584	1470	1515	1482
Mass density [kg/m${}^3$]	1150	1099	937	1040	1000

,

Table 2. Thermodynamic parameters corresponding to various parts of breast tissues used for simulation of heat deposition, diffusion, and the subsequent calculation of lesion area [36].

Characteristics	Skin	Vessel	Fat	Gland	Water
Thermal conductivity [W/(m·K)]	0.37	0.46	0.21	0.33	0.6
Specific heat [J/(kg·K)]	3391	3306	2348	2960	4178

and

Table 3. Optical parameters used for simulation of light propagation for initial pressure distribution via photoacoustic effect [35].

Characteristics	Skin	Vessel	Fat	Gland	Water
Scattering coefficient [cm${}^{-1}$]	500	179	159	133	0
Absorption coefficient [cm${}^{-1}$]	0.08	9	0.05	0.04	0
Scattering anisotropy	0.99	0.975	0.95	0.95	0.99
Refractive index	1.4	1.38	1.4	1.4	1.33

. The region surrounding the breast (the black region in Figure 1) corresponds to water that serves as an acoustic coupling medium.

2.2. Ultrasound Propagation

Each transducer in HIFU is a source of sinusoidally propagating ultrasonic waves (at higher intensity) to the target tissue medium via a coupling medium (water in our case). Mathematically, we can write the wave equation governing the propagation of acoustic waves as [37]:

$${\nabla }^{2}p(r,t)-\frac{1}{C_s^2}\frac{{\partial }^{2}p(r,t)}{{\partial }^{2}t}=0,$$

(1)

where p and $C_s$ are the pressure and ultrasound speed, respectively. Acoustic waves are sensitive to reflection, transmission, and refraction at an acoustically varying tissue–tissue interface. Because of this mismatch in acoustic parameters, a fraction of the acoustic wave energy is reflected at every tissue interface, while the remaining energy is concentrated at the HIFU focal point inside the tissues. In this way, ultrasonic waves undergo attenuation while propagating through a mechanical medium. In contrast to the case of fluids, the acoustic wave equation for tissue media obeying a power law for pressure attenuation is governed by the causal wave equation (which differs from Equation (1) above). It is fundamentally derived from the fractional constitutive equation between stress and strain [38]. The attenuation coefficient ($\alpha$), which is defined as $\alpha ={\alpha }_a+{\alpha }_s$ (where ${\alpha }_a$ is the absorption coefficient and ${\alpha }_s$ is the scattering coefficient), provides a measure of the loss of acoustic energy in a medium as it propagates. Following a power law (as is the case of tissues), the attenuation coefficient in terms of ultrasonic frequency is expressed as [39]:

$$\alpha =a{f}^b,$$

(2)

where a and b are constants, and f is the ultrasound frequency. Presuming attenuation is constant over a distance (x), the pressure amplitude decays exponentially with the distance traversed by the wave as a result of these attenuations, as expressed by [39]:

$$P_2=P_1{e}^{-\alpha x},$$

(3)

where $P_1$ and $P_2$ are acoustic pressures at two points inside the tissue separated by a distance of (x) along the direction of wave propagation. As a result, the acoustic (harmonic plane) wave solution is altered [39]:

$$P(x,t)=U_0{e}^{-\alpha x}{e}^{i\omega (t-\frac{x}{C_s})},$$

(4)

where $U_0$ is the pressure at $t=0$, and $x=0$.

For the simulation corresponding to HIFU, we used k-Wave [31], which is an open-source MATLAB toolbox, for the time-domain simulation of acoustic wave fields. A snapshot of the simulation of acoustic wave propagation overlayed on the cross-section of an extremely dense-type breast is shown in Figure 2a, along with the distribution map of mass density (Figure 2b) and the acoustic attenuation coefficient (Figure 2c). The frequency was fixed at 1 MHz throughout the simulation studies, and the HIFU was focused approximately at the same distance from the skin tissue in all three breast samples to maintain uniformity. In Figure 2 for the case of an extremely dense-type breast, the simulation area was 92.25 × 112.8 mm${}^2$ with a pixel size of 0.2 mm × 0.2 mm. A HIFU arc length of 51 mm with a radius of 40 mm was simulated by a number of point pressure sources (∼339 with a pressure amplitude of 0.7 MPa) arranged uniformly on an arc.

2.3. Heat Deposition and Diffusion

The lost energy from ultrasound attenuation gives rise to local heating within the soft tissue [39]. The rate of heat deposition ($(Q)$) is [40]:

$$Q=\frac{\alpha P^2}{\rho C_s}$$

(5)

where $\alpha$ is Nepers per length, and P is pressure amplitude. The deposited heat is not localised with time because of heat diffusion through thermal conductivity (as well as blood perfusion). The transient heat transfer equation provides a general equation controlling thermal dynamics within tissues, and it can be written as [39]:

$$\rho C\frac{\partial T}{\partial t}=k{\nabla }^{2}T-A(r,t)+Q,$$

(6)

where k, C, and T are thermal conductivity, specific heat, and absolute temperature, respectively. This is basically an energy conservation statement, where $\rho C\frac{\partial T}{\partial t}$ is the net heat deposition rate per unit volume; $k{\nabla }^{2}T$ is the thermal diffusion term; $A(r,t)$ accounts for vasodilation and perfusion, which serves as a sink (cooling source); and Q is the heating source (from HIFU). Upon neglecting the cooling term in Equation (6) and using Equation (5), we obtain:

$$\rho C\frac{\partial T}{\partial t}=k{\nabla }^{2}T+\frac{\alpha P^2}{\rho C_s},$$

(7)

In order to non-invasively and non-destructively ascertain the heat deposition and temperature rise from HIFU, we used the above equation for numerical simulation. In this way, we adopted Equation (5) for numerical simulation of the rate of heat deposition using the k-Wave toolbox. The thermodynamic parameters adopted for this part of the simulation are shown in Table 2.

2.4. Thermal Dose and Lesion Formation

In order to access thermal tissue damage during heat deposition, one must determine the thermal dose received by the tissue. A unit called cumulative equivalent minutes ($CEM43$) is commonly used to convert various time–temperature exposure profiles in experiments into an equivalent exposure time expressed in minutes with a reference temperature of 43 °C [41]. Mathematically, it is expressed as:

$$CEM43 °\mathrm{C}={\int }_0^t{R}^{43-T}dt,$$

(8)

where $R=1/4$ for $T<43$, $R=1/2$ for $T>43$, and T is the temperature in °C.

In a wide range of tissues, regions with $CEM43>240$ generally result in permanent ablation. We consider tissue regions with $CEM43>240$ to be the area of the lesion in our current k-wave simulation studies. Due to heat diffusion in a medium, the lesion typically progresses from the focal point of the HIFU on-time. In the k-wave simulation, the lesion area was considered an area of grid points where $CEM43$ equalled or exceeded 240.

2.5. Initial Pressure Rise in PA Imaging

The phenomenon of the generation of an acoustic signal resulting from the tissue’s rapid heating brought out by the transient illumination and absorption of a ∼nanosecond pulse of electromagnetic radiation is referred to as the photoacoustic (PA) effect. The lesion provides excellent contrast in the reconstructed PA image due to the contrast in the distribution of the lesion’s optical absorption coefficient with respect to the surrounding normal tissue. The rapid heating caused by transient electromagnetic wave irradiation can be seen as a perturbation of the thermodynamic equilibrium [28], since the tissue system is completely characterized by pressure ($(P)$), volume ($(V)$), and temperature ($(T)$). Therefore, we can consider the tissue system to be a hydrostatic thermodynamic system. The initial pressure wave (initial photoacoustic signal, $(P_0)$) is related to the local temperature rise ($\Delta T$) at the lesion site that is caused by a short laser pulse and can be expressed as [28]:

$$P_0\approx \frac{\beta }{K_T}\Delta T,$$

(9)

where $\beta$ is the volume expansibility, and $K_T$ is the isothermal compressibility. Using the thermodynamic heat transfer relation for the absorbed heat density ($\Delta Q_{dens}=\rho c_v\Delta T$) and assuming that the optical energy absorbed by the tissue is completely converted into heat energy so that $\Delta Q_{dens}={\mu }_a\varphi$, we have

$$P_0=\frac{\beta }{K_T}\frac{1}{\rho c_v}{\alpha }_a\varphi ,$$

(10)

where $c_v$ and $\varphi$ are specific heat capacity at constant volume and optical fluence, respectively [28], and $\Gamma =\frac{\beta }{K_T}\frac{1}{\rho c_v}$ is a tissue-specific constant at a particular temperature known as the Grueneisen parameter. The tissue medium then acquires this initial pressure distribution, which is represented by $P_0(r)(=\Gamma {\alpha }_a\varphi )$.

We employed ValoMC [32], a Monte Sarlo software (MATLAB toolbox), to simulate the transport of visible and near-infrared-range photons in biological tissue to simulate the photoacoustic effect and induce an initial pressure rise distribution. All the optical parameters used are shown in Table 3. The Gruneisen parameter was kept constant at 0.02. A Gaussian beam with a diameter of 30 mm and a standard deviation of $\sigma =0.1$ was employed. The beam direction is shown in Figure 3a, as indicated by the red arrows from the bottom. The absorbed optical energy density was calculated by taking the product of optical fluence ($\varphi$) and the optical absorption coefficient (${\alpha }_a$) at each point. Finally, the initial pressure rise was calculated by taking the product of the Gruneisen parameter and the calculated optical energy density. Note the exceptional visibility of the blood vessels in the initial pressure rise distribution map (Figure 3b) resulting from a very high relative value (hundreds of orders higher than the other included tissues) of the absorption coefficient (${\alpha }_a$) for the vessels.

2.6. Reconstruction

An assortment of acoustic sensors acquires the initial PA waves ($P_0$) upon exiting the tissue boundary. These boundary-collected data are used to reconstruct PA images using a variety of reconstruction methods. In the simulation study, a simple delay-and-sum (DAS)-based beamforming algorithm was used for the reconstruction of internal pressure distribution ($P_0$) resulting from transient illumination. As shown in Figure 4a, $P_0$ was acquired by a linear array transducer shown by a white line at the bottom. Figure 4a is the same as Figure 3b, except that the color range was squeezed to bring out more details. Figure 4b shows the corresponding reconstructed image. The sampling frequency was 33.33 MHz.

3. Results

The simulation results of HIFU focused on three different types of breast (extremely dense, heterogeneously dense, and scattered fibroglandular) are shown in Figure 5. Our aim here is to determine how precisely the ultrasound waves from the HIFU are focused at its focal point. To this end, all the parameters of HIFU, such as arc length, radius, etc., were kept identical for all three breast types (except for HIFU position and tilt so as to have the freedom to align with varying shapes or sizes of breasts). The distance between the HIFU arc and the skin tissue was also kept the same for all breast types. The HIFU parameters are the same as those described in Section 2.2. The first column (Figure 5a–j) corresponds to the extremely dense-type breast, the middle column (Figure 5b–k) corresponds to heterogeneously dense-type breast, and the last column (Figure 5c–l) corresponds to scattered fibroglandular-type breast. The top row (Figure 5a–c) shows the distribution of HIFU pressure wave amplitude overlayed on the respective breast type. As can be seen, it is maximum for scattered-type breast at the focal point of HIFU, where the pressure waves are sharply focused. The second row (Figure 5d–f) shows the distribution of the heat deposition rate per unit area. Due to the relatively high presence of fibroglandular tissue in the focal zone in the case of heterogeneously dense breast, the heat deposition rate is maximum among the three breast types (see Section 4 related to Figure 6). The third row (Figure 5g–i) shows the lesion map overlayed on the respective breast type. It is interesting to note that despite the varying structures of tissues within, the distance from the centre of lesions to the HIFU arc came out to be the same for all three breast samples. The bottom row (Figure 5j–l) shows the initial pressure ($P_0$) distribution using ValoMC for the induced lesion. For this, the blood vessels were suppressed, as they have optical absorption coefficients hundreds of orders of magnitude higher than that of the surrounding breast tissue, resulting in high pressure at those sites.

Figure 6 shows the distribution of the heat deposition rate in the extremely dense-type breast. Upon close analysis of the heat deposition distribution and the overlayed breast image, we can see that the rate peaks in the fibroglandular tissue if it happens to be in the vicinity of the focal zone of the HIFU.

The variation of lesion area (area with $CEM43>240$) with HIFU on-time for the three different types of breast is shown in Figure 7. Again, the HIFU parameters were kept the same, as described in Section 2.2 while increasing the HIFU on-time. Figure 8 shows a graph demonstrating the temperature increase exactly at the focal point of HIFU within the breast with HIFU on-time (5 s) and off-time (60 s). The temperature was measured at time steps of 0.1 s. This was done for all three samples. The peak temperature for the scattered-type breast came out to be the maximum among all, as expected, because of the obtainable sharp focusing (thus, more HIFU fluency) at the focal point compared to the other two samples.

4. Discussion

HIFU can be a viable, safe, effective, and non-invasive therapy for the treatment of localized breast cancer without any serious negative side effects [42]. Here, a 2D simulation for light- and acoustic-field propagation was carried out on three sample types of contrast-enhanced MR breast—scattered, heterogeneously dense, and extremely dense type (Figure 1). After incorporating the distribution of various required physical properties (Figure 2), the simulation for the HIFU field was carried out. The resulting thermal dynamics with lesions generation are discussed in Figure 9 and Figure 10. The initial pressure distribution determined via photoacoustic effect was also simulated (Figure 3). The results show that although breasts are made up of soft tissues, the complex structures within them do not show the propagation of acoustic waves similar to those of other soft (homogeneous) tissues. Our numerical simulation showed that scattered-type breasts can benefit most from HIFU therapy.

From Figure 5d–f, we can infer that the denser the structure of fibroglandular tissue within the breast, the less sharply the ultrasound is focused at the focal point of the HIFU. Ideally, the lesion should have been elliptical in shape [43]. The deviation from the ideal behaviour is due to attenuation and scattering of focused ultrasound as it propagates through the breast, where it encounters fat and glands with different acoustic properties, which are distributed non-uniformly inside the breast. This imposes a defocusing challenge for HIFU therapy. Thus, HIFU is less precise for extremely dense breasts, as can be seen in Figure 5d. Our results are in agreement with those reported in [44], in which researchers studied a way to reduce defocusing by focus control based on the time-reversal method.

The initial pressure distribution shown in Figure 5j–l was carried out using the ValoMC MATLAB toolbox [32] to simulate photon transport within the tissue. The direction of the Gaussian beam is evident from the bottom, where intensity is high. The blood vessels were masked by the properties of fat tissue for this particular case to explicitly study pressure distribution at the lesions. The lesion areas were given a value of optical absorption coefficient 0.5 cm${}^{-1}$ higher than that of fat [45]. As expected, the magnitude of pressure falls as we move along the propagation path of the light beam. The pressure magnitude at the focal point (lesion centres) for scattered fibroglandular breasts is higher than that of heterogeneous and extremely dense-type breasts. However, because of an extended lesion area due to defocusing in extremely dense breasts, the tip of the lesion comes closer to the direction of the light source, which results in a higher magnitude of pressure in that part. Thus, tumours or lesions can be detected relatively easily in scattered fibroglandular-type breasts.

The map of heat deposition rate per unit area (Figure 6) hints at yet another problem we might encounter while performing HIFU therapy. Parts of the fibroglandular tissue that are in the vicinity of the focal point of HIFU where the heating is to be performed acquire a locally higher rate of heat deposition. The reason can be traced to the acoustic attenuation coefficient relative to the surrounding fat tissue. Depending on the size and location of the tumour, this can pose an issue of potentially harming healthy glands in the vicinity of the HIFU focal point. This may not pose much concern in the case of scattered fibroglandular-type breasts, where the HIFU focal point can be sharply focused with little intervening fibroglandular tissue structure. So far, we have not found any source, either simulation or experimental, to verify this specific result. This can be a subject for future studies.

Finally, an investigation into the development of lesion area and temperature with HIFU on-time was carried out. This study demonstrates (see Figure 7) that we can assume that the lesion area grows fairly linearly over a time scale of a few seconds depending on the HIFU operational parameters. This is in agreement with lesion imaging experiments carried out in soft tissues and reported in Ref. [46]. For the given HIFU parameters and on-time, the scattered-type breast is observed to have a lower rate of increase in lesion area compared to extremely or heterogeneously dense breasts. This may be because of the presence of fibroglandular tissues in the vicinity of the focal point of the HIFU in the latter two samples, where the temperature peaks sharply (thereby essentially spreading the lesion centres), as discussed in the preceding paragraph. The relationship between lesion size and HIFU exposure time is not strictly linear for high-intensity and long time scales. While longer exposure times generally result in larger lesions, the relationship is complex and depends on several factors, including the intensity and frequency of the ultrasound waves, as well as the tissue type and properties. In general, there is a threshold exposure time below which HIFU treatment may not result in a visible lesion. Beyond this threshold, lesion size tends to increase with exposure time. However, there may be diminishing returns beyond a certain exposure time, as tissue may become more resistant to HIFU-induced thermal damage. Additionally, lesion size may not increase linearly with exposure time in all cases. Factors such as tissue heterogeneity, blood flow, and tissue perfusion can affect the ability of HIFU to induce tissue damage and may result in non-linear relationships between exposure time and lesion size. For a theoretical prediction of the size of a thermal lesion in soft tissue during HIFU treatment, see [47]. There have also been integrations of external controllers to achieve preferred and linearly varying predictable lesion sizes. For an experimental work in which the temperature was controlled to achieve the same, see [48]. Figure 8 shows the temperature increment and decrement at the focal point with HIFU on-time (5 s) and off-time (1 min). These trends are also in line with experimental (phantom) results (see Ref. [49]) and numerical simulation [50] for homogeneous soft tissues. The expected temperature at the focal point of scattered-type breast at any time is observed to be higher than for the other two types due to sharp focusing. Due to limited availability of hardware resources in our lab, experimental validation studies were not performed in the present work. We reserve this as an aim for our future research studies.

This shows that acoustic inhomogeneity plays a central role in the effectiveness of HIFU therapy in breasts. An optimal arrangement with respect to an individual’s structural composition of breast and identification of suitable HIFU operational parameters via a simulation such as this will greatly aid in achieving the goal. From the perspective of practical (realistic) applications in the present study, one can obtain a 3D map of the patient’s breast from a quality imaging modality (such as MRI) and then load it in a simulation software to virtually simulate HIFU therapy as in the present study. Such simulations, being specific to a patient’s breast structure and composition, can be of great impact in deciding the HIFU operational parameters and their alignment. In future works, we plan to incorporate the reconstruction of initial pressure distribution for photoacoustic imaging to monitor HIFU therapy in real time. The reconstruction of such lesions inside the breasts, along with intervening acoustic inhomogeneity, would require an advanced beamforming algorithm to have some correlation with lesion shape and size instead of a simple DAS beamforming algorithm (see Figure 4b for the sake of completion). The inclusion of machine learning will a great aid on this front.

The study presented in this article is not limited and can be further improved by incorporating several other aspects. In conventional simulations, the tissue parameters are kept constant for predictions, as in the present study. However, accounting for the change in the acoustic absorption coefficient due to changing temperature at the focal point leads to higher focal zone temperature variation and larger lesion sizes [43]. On the other hand, the incorporation of blood perfusion will serve as a sink for thermal energy at the focal zone, thereby reducing the focal zone temperature change and therefore the lesion size. Another aspect to take into account, which is more like a challenge, is the possibility of acoustic cavitation in the focal zone. Vapour bubbles could be induced if the temperatures at the focal zone approach the boiling point. This leads to enhanced heating and contributes significantly to induced hypothermia [51]. It is interesting to note that, as we have shown here, temperature may have a sharp peak in fibroglandular tissue if it lies in the vicinity of the focal point of HIFU; consequently, these sites can be more susceptible to bubble formation. There is also an effect on the magnitude of initial pressure distribution (from the photoacoustic effect) over the lesion created by HIFU. Both the optical absorption coefficient and reduced scattering coefficient increase in the lesion area compared to native tissue [45,52]. Better segmentation of breast tissue structures with more accurate values of their physical properties will obviously enable an overall improvement.

5. Conclusions

In this study, we simulated HIFU therapy by employing a MATLAB-based k-Wave toolbox and generated initial pressure distribution from the photoacoustic effect using the ValoMC toolbox. Both the open-source numerical tools—k-Wave and the ValoMC—were integrated together in MATLAB for the processing of three different numerical breast phantoms—extremely dense, heterogeneous, and scattered fibroglandular-type breast—along with the details of theoretical (analytical) formulation. We discussed the implications of the results and the associated challenges for HIFU therapy. We concluded that scattered fibroglandular-type breasts are well-suited for HIFU therapy and even for imaging of lesions or tumours due to their low acoustic inhomogeneity. Further improvements in the simulation can be achieved by incorporating dynamic tissue property variables with changing temperature and performing the simulation in 3D. Employing a better segmented numerical breast phantom, along with appropriate sources for acoustic and optical property values, will increase the reliability of the simulation. In the future, we plan to perform experimental validation studies and incorporate machine learning for the reconstruction of initial pressure distribution for photoacoustic imaging to monitor breast HIFU therapy in real-time.