Original paper
Multi-axis dose accumulation of noninvasive image-guided breast brachytherapy through biomechanical modeling of tissue deformation using the finite element method

Purpose

Brachytherapy (BT) is an established treatment modality that places the radiation source within or in close contact with the patient [1]. Through the positioning of needles or applicators, it is possible to deliver dose distri­butions with BT that are more conformal than possible with external-beam radiotherapy (EBRT) using X rays, protons, or carbon nuclei [2]. Modern treatment planning techniques utilize images from computed tomography (CT), magnetic resonance imaging (MRI), and/or positron emission tomography (PET) for performing image-guided dose calculations. For EBRT, regions of interest such as healthy tissues or the target are contoured from these images, while CT values (i.e., Hounsfield units) are used quantitatively for dose calculations. Historically, BT has used images only for contour delineation, with dose calculations based on the superposition of single-source dose distributions in a liquid water environment [3,4]. Recently, images have also been used for BT dose calculations [5,6], where dose calculation errors without accounting for tissue material heterogeneities have been shown to range from < 5% to over a factor of ten [7].
With improved accuracy for BT dose calculations now in the spotlight, the clinical process of image utilization is being reevaluated [8]. Image registration is the key step in image-guided radiotherapy. Yet our approach to image utilization has been simplistic for the past two decades, assuming that the patient’s anatomy does not change between the simulation and treatment, over the course of radiation therapy, and is a rigid body [7]. When attempting image fusion with multiple imaging modalities, such as MRI or PET to a reference CT dataset, it quickly becomes clear that the patient is setup in a different pose and that simplistic approaches based on image translation or rotation will not permit overlap of all relevant parts of the patient’s anatomy, even when restricted to the most local tissues [8,9].
Deformable image registration (DIR) is a computational method for comparing multiple image sets for the same subject, yet having different geometries [10,11]. Simply, volumetric information from one imaging modality is warped in a robust manner to mimic the anatomical geometry of the principal imaging modality. In addition to providing the ability to integrate data from different imaging modalities (i.e., CT, MRI, and PET), DIR is used in radiotherapy treatment planning for motion assessment (e.g., 4D CT scans) and also for dose accumulation across multiple treatment fractions or treatment modalities [12]. In this latter manner, there is promise for adaptive radiotherapy during a given treatment course where imaging data from prior treatment fractions are used to guide subsequent treatment fractions, instead of delivering the same treatment plan from start to finish [13,14].
In the context of BT, DIR has been used for co-registration of images from CT and MRI for improved target delineation [15,16], for dose summation of separate treatment fractions [17,18], and for dose summation of EBRT and BT treatments [19]. This can be especially important for the accurate dose calculation of high-energy photon-emitting BT sources [20]. Several DIR algorithms are commercially available. These include MIM SymphonyTM by MIM Software, Inc. (Cleveland, OH, USA), VelocityTM by Varian Medical Systems, Inc. (Palo Alto, CA, USA), RTxTM by Mirada Medical (Oxford, UK), and IMSimQATM by Oncology Systems, Ltd. (Shropshire, UK). There are also numerous other systems that are developed through academic research [21].
Deformable image registration algorithms may be rank­ed based on their utilization of biophysical properties, starting with a simple dataset translation, translation + rotation, affine (e.g., translation + rotation + shear), spline techniques (e.g., thin-plate spline for global adjustments and b-spline for local adjustments), and biomechanical models using material properties (e.g., tissue compressibility, stiffness, boundary effects). This latter type of DIR algorithm aims to go beyond the morphological limitations and arbitrary imaging choices that govern the other DIR approaches.
Among several treatment modalities for partial breast irradiation, noninvasive image-guided breast brachytherapy (NIBB) delivers conformal treatments using BT applicators aligned to the target using mammographic imaging [22,23]. First-in-class NIBB applicators, AccuBoost® [24-26], are available from Advanced Radiation Therapy, LLC (Tyngsboro, MA, USA), and have been used for treating the tumor bed as a radiation boost for early-stage breast cancer as well as for accelerated partial breast irradiation (APBI) monotherapy. While initial clinical results have been favorable [27,28], a practical clinical challenge has been to depict summed dose from each set of parallel-opposed BT beams [29]. The collimated dose distribution from individual applicators can be incorporated into a conventional BT treatment planning system [30]. However, dose summation from treatments on (or near) orthogonal axes is not feasible given the extent of orthogonal breast compressions intrinsic to the modality. Further, mammography is used for high-contrast, high quality images to locate the tumor bed and to position the applicator over the target. However, it is not possible at this time to gather 3D volumetric imaging such as from CT or MRI during the time of applicator positioning and treatment delivery [31]. Consequently, there is a need for a site-specific (i.e., breast) DIR model that does not rely on imaging. The purpose of the current study is to develop a biomechanical model of the breast that can account for orthogonal compressions and dose summation as delivered using NIBB.

Material and methods

A 2D section of breast tissue was modeled as a finite element mesh with biomechanical properties obtained from the literature. A finite element analysis (FEA) was applied to the model as it stepped through varying levels of compression. In order to utilize this model for a 3D geometry, it was assumed that the breast was an infinite right cylinder. This approach assumed there were no interactions between successive sections and allowed 2D sections to be stacked sequentially an arbitrary number of times to build up the 3D model of the breast. Boundary conditions at the chest wall and those due to skin effects at the anterior aspect of the breast were ignored. Dose distributions from the collimated beams of the standard round and the skin-dose optimized (SDO) applicators have cylindrical symmetry. Results from each Monte Carlo (MC) simulation of the dose distributions are provided as a 2D section from a single applicator source. To model the 3D dose distribution for a set of opposing applicators, symmetry from the MC data was utilized to obtain the full 3D dose distribution. This dose was then applied to the biomechanical model for a given compression in both the cranial-caudal (CC) and medial-lateral (ML) directions to determine the dose delivered to each location within the tissue. The dose distribution throughout the breast was then available for determining a number of treatment planning parameters.

Biomechanical model

Biomechanical modeling of the breast deformation during mentioned procedures (NIBB, AccuBoost, APBI monotherapy or mammography) is non-trivial. Complex anatomy, uncertainty in biomechanical material properties, and inhomogeneity of different tissue types (fat, fibro- glandular, lesion, ligament, and skin) within the breast contribute to the challenges of an accurate model. Additionally, such models are expected to simulate real-life response and behavior of the breast in procedures such as biopsy and mammography that involve large deformations of the organ. The presence of very soft tissue, such as fat, poses another difficulty as classical FEA modeling has mathematical limitations in cases of extreme deformations.
Research efforts cover a wide range of approaches to the mentioned challenges, which include assumptions of homogenous material [32-34], incompressible vs. compressible tissue [35,36], more detailed modeling of the internal structure of the breast using information from MR scans [37], non-linear material properties [38,39], and hybrid methods that combine Lagrangian mechanics with the Eulerian method [40].
An increase in complexity of the modeling approach, which may result in higher accuracy of the deformation simulation, has a direct impact on the computational requirements of the analysis. Therefore, one must choose a biomechanical model that balances predictability and accuracy requirements against time and setup complexity.
Given the high deformability of the breast and the large axis deformation imposed during clinical NIBB procedures, a complete form of finite elasticity theory should be considered as the starting point in the modeling stage. However, some assumptions are needed in order to implement a practical solution of the theory for this project. First, it was assumed that the material is isotropic hence its mechanical properties are identical in three dimensions. Additionally, it was assumed that the breast is under planar stress. This assumption is based on the idea that the breast tissue is less constrained to move orthogonally with respect to chest wall, thus the stress in that direction can be assumed to be negligible, which is in line with the premises of a plane-stress loading. Experimental results have also been in favor of plane-stress approach as opposed to plane-strain ones [41]. Considering the above issues, the partial differential equation (PDE) describing the displacements can be
simplified as: E U V E U V ––––––– (– –––––– + v –––––––– )+ ––– ––––––––– ( –––––– + –––––––– )= 0 Eq. (1) 1–vxxy2 (1+v) yxy E V U E V U ––––––– (––––––– + v –––––––– )+ ––– ––––––––– ( –––––– + –––––––– )= 0 Eq. (2) 1–vyxy2 (1+v) xxy
where U and V are displacements, E is Young’s modulus of elasticity, and v is Poisson’s ratio. Eq. (1) and Eq. (2) essentially describe the relation between applied force and resulting deformation in an elastic medium. This approach to modeling biomechanical behavior of tissue as a stress-strain problem has previously been used to account for tumor and tissue deformation. For example, Li et al. employed spatially-varying E values in a hybrid FEA and DIR approach to account for respiratory induced deformation to improve the accuracy of calculated dose distributions in thoracic radiation therapy [42].
Using the Galerkin method of weighted residuals as the FEA strategy and assuming constant strain elements [37], a model was developed in MATLABTM (version 2012a, The MathWorks, Inc., Natick, MA, USA) to solve Eq. (1) and Eq. (2) over the domain. Using FEA, the problem domain was discretized into small regions (consisting of nodes and elements), over which the PDE is solved. This essentially transforms the above equations into a system of linear equations in the form of [F–] = [K] · [U–], where F– is nodal force vector and U– is the nodal displacement. The matrix K contains information about material properties of the elements and is called the stiffness matrix, which is a function of E and v. The former parameter dominates tissue movement within the plane of analysis while the latter dominates tissue flow out of the plane of analysis.
The domain of the problem (breast) was modeled as a circle with a 15 cm diameter, which was subsequently partitioned using Triangle software (Berkeley, CA, USA) [43] to produce a 0.1 cm thick planar mesh with a fixed number of about 1,300 triangular elements and 700 nodes (Fig. 1A). These planar meshes were adjacently stacked to subtend the 3D volume.
As previously described, a range of biomechanical properties of tissue has been reported. Material properties of the breast are mainly calculated using experimental approaches such as elastography [44] and diffuse wave spectroscopy [45]. Van Houten et al. [46] have compiled a list for E values for normal fat and normal glandular tissue. They used values reported by others as well as providing results from their own elastography study. As reported by Van Houten et al. [46], the average value of E = 30 kPa for healthy tissue was selected for use throughout the entire volume of the current study. As illustrated in Figure 1A, the model used in the current study is assumed to have homogenous material properties.
Compressibility of the breast tissue is another important parameter in modeling step, which is dictated by v. Values 0.45  v  0.49 are representative of an incompressible material. However, since fatty tissue is not completely incompressible, a value of v = 0.30 was selected for the current study [46].
Since the deformations that are applied to breast using compression plates are substantial, the displacements of the elements can lead to unpredictable results and computational divergence. To model this behavior, we divided the large motion of the plates into smaller steps of 0.01 cm. Within each step, since the deformations fall into the confinements of linear elasticity, we can safely use the methods and equations described previously to calculate the strain using Cauchy’s strain tensor formula [47]. The simplicity and relative ease of implementation of this method results in rather fast computations.
Additionally at each step, the boundary conditions of the elasticity PDE require updating since a larger area of the breast would come into contact with compression plates. This is accomplished by arresting all the degrees of freedom of the perimeter nodes in contact with plates. The increasing number of arrested nodes as the algorithm runs through each step essentially provides the nodal forces that populate F–. The data obtained from solution of PDEs as outlined above represents the nodal displacements of the vertices within the mesh in every step of the breast deformation until it reaches its final compressed state.

Mapping dose accumulation

For this analysis, the breast was approximated as an infinite right cylinder, separated into 0.1 cm-thick circular sections from the FEA described above. The FEA output coordinates were gathered for the desired compression case to establish the dose delivered to each unit of tissue. A 3D MC dataset was queried to determine the dose at the coordinates of each node in the compressed state. That dose was then tied to the unique nodal identifier and held as the nodal coordinates reverted to the uncompressed case. This process with applicator beams oriented along the CC direction was then repeated with a 90° rotation along the X-axis to model the subsequent ML half-fraction.
As described in Rivard et al. [25], the AccuBoost applicators, round or SDO [26], were geometrically modeled using radiation transport simulations to obtain dose distributions in breast tissue as a function of applicator diameter and compression thickness. The resultant dose distribution from each applicator is cylindrically symmetric about the applicator central axis due to the applicator design. The grid size was 0.1 cm in both the depth and radial directions from the applicator central axis, and subtended a height of 12 cm and width of 15 cm.
A single axis, dual-vector, dose delivery consists of two opposing coaxial applicators on either side of the breast. In clinical applications, the compressed breast thickness ranges from 3 cm to 8 cm. The range used in this analysis was consistent with the APBI guidance. To simulate the dose delivered by the second applicator, the array was truncated at the desired breast thickness and mirrored about the midplane of the compression plates. The individual dose fields were then superposed upon one another to determine the dose distribution (see Fig. 1B).
2D data matrices were created in this manner for all possible combinations of separation, applicator size, and applicator type. Identical applicators were used in each orthogonal compliment dual-vector dose delivery. 2D dose distributions were created in the CC orientation and rotated 90° for ML distributions. Clinically, 1 Gy is prescribed to the center of the breast per axis. A dual axis treatment of 1 Gy delivered CC (Fig. 1B) and 1 Gy delivered ML (Fig. 1D) was used. Since the 1 Gy from the two beams do not coincide for non-centered targets, the beams must be normalized separately rather than scaling to the center of the summed result. At this point, all matrices were normalized to their central value and scaled to 50%. This facilitated a normalization, where 100% corresponds to 2 Gy for properly handling of centered as well as non-centered targets.
The algorithm was used on all the nodes until the 3D field delivered to the breast in the compressed state was filled. The geometry of the round applicators allows for parts of the matrix to be mirrored, avoiding redundant calculations of symmetric volumes. Using a primary plane radius of 15 cm, the resulting matrix represented a 30 cm diameter cylindrical dose distribution with varying height based on compression separation, though much of the dose distribution outside a 15 cm diameter is beyond the simulated breast volume.
Each XY plane of the 3D dose distribution was applied to the 2D FEA model. The algorithm processed each plane in the dose distribution sequentially from +Z to –Z. In applying the FEA model, the XY coordinates of each node were overlaid on the dose distribution plane. Figure 1E shows the resultant dose distribution plane and nodal coordinates using 5 cm round applicators with prior breast compressions of 5 cm in both directions. Bilinear interpolation of the dose field was used to find the dose delivered at each set of nodal coordinates in the compressed state. The figure illustrates the beam arrangement for a centralized target.
The dose matrix for the uncompressed breast cylinder was created using the dose delivered to each node along with its uncompressed coordinates. The MATLAB function TriScatteredInterp was used to interpolate dose intensities between uncompressed nodes into 1 mm cubic gradations. Cranial-caudal and ML axis were simulated separately and added together to create the full ortho­gonal beam fraction. Given the size of the FEA model, the uncompressed matrix was 15 cm by 15 cm by 30 cm along the X, Y, and Z respectively. This matrix represents the dose delivered to a 7.5 cm radius right circular cylinder breast model. As measured in the Z-axis, the height of the cylinder is arbitrary, so a value of 30 cm was used to encompass the entire 30 cm-wide dose distribution and to assure the entire compressed breast volume remained within the model coordinates.
With the relative magnitude of the dose delivered to every node in the 3D breast model being determined, dose value histograms (DVH) and dose contour plots were developed for each set of parameters. The DVH plots depict the percentage of a spherical target volume that has a dose value above a given value. This entails sampling a set of internal nodes within a specified radius of the target location and quantifying the dose relative to the 2 Gy prescription dose. Another clinically useful parameter is the dose applied to the skin in the vicinity of the applicators relative to the center dose. The external nodes directly under each applicator were identified and their dose sampled. The average dose of the sampled skin was then normalized to the 2 Gy prescription and tabulated with the same method as used by Yang and Rivard [26].

Model deformation and the affect of compression on target size

Table 1 shows the affect of compression thickness on target size in the planes perpendicular to the axes of compression. This change in target size with compression needs to be accounted for appropriate applicator selection. Targets were also simulated with varying offsets from the center of the breast model. For dosimetric coverage of the target, the applicator field was consequently shifted. These offsets were in the positive X direction (along the ML direction) and perturbed the initially spherical shape of the target as shown in Table 1.
Based on info in Table 1, these constraints were used to guide applicator choice for subsequent evaluations of dose distributions and DVHs, where plots were generated using the uncompressed distribution matrices. Additional plots of the dose distributions on the skin were obtained. Three spherical targets were modeled with diameters of 2.5 cm, 3.5 cm, and 5 cm. To calculate the DVH, dose values that fell inside the target sphere were isolated into an array, where each value represented a (0.1 cm³) volume. The intensity and volume was used to plot the DVH curves. The dose was not normalized prior to plotting. Consequently, the 100% isodose line represents 2 Gy, where 1 Gy was prescribed to the center of each compression. Targets that were not centrally positioned may result in DVHs, where simulated tissue received more dose than at the center.

Benchmarking from original Monte Carlo

Towards supporting the methods presented in the current study, comparisons were made with the published MC data. The MC results of the skin-to-center dose ratio (SCD) for the 5, 6, 7, and 8 cm diameter round and SDO applicators [25,26] for compressions of 3 cm to 8 cm were compared to SCD results from the current study preceding summation of the orthogonal set of parallel-opposed beam, which did not include such “scatter” contributions. For the 24 conditions examined for the round and SDO applicators, the SCD ratios from the current study to the published results were 0.993 ± 0.007 (min. 0.979, max. 1.004) and 1.004 ± 0.006 (min. 0.991, max. 1.017) for the round and SDO applicators, respectively. Considering the effect of volume averaging on the 0.1 cm-thick voxels and the assumed values of v and E, these dosimetric differences of no more than 2% were considered acceptable.

Results

In the following sections, the FEA model is systematically evaluated in terms of target coverage and skin dose in relationship to several variables over a clinically relevant range.

Effect of round applicator size on target coverage and skin dose

Applicator size had only a small effect on target coverage until the applicator size was less than the target size at compression. In all cases where an applicator larger or equal to the target size at compression was used, resulted in > 90% of the target covered by > 90% of the prescription dose. Figure 2 shows representative DVHs for 3.5 cm and 5.0 cm diameter uncompressed targets irradiated at 6 cm compression with varying diameter round applicators. From Figure 2A, it is evident that the smaller target is covered sufficiently by all the possible applicator sizes. Figure 2B demonstrates significant underdosing that would occur when irradiating a 5 cm diameter target (6 cm compression) using a 5 cm diameter round applicator. This situation improves as the applicator size increases. However, care must be taken to concurrently minimize irradiation of healthy tissues. This process of assessing target coverage is representative of current clinical practice.
Examining target positions with offsets of 0, 2.375, and 4.75 cm from the central positions, the DVHs for covering 3.5 cm and 5 cm diameter uncompressed targets with 6 cm and 8 cm diameter round applicators, respectively, are shown in Figure 2C and Figure 2D, respectively. In all cases, the dose coverage became less uniform as the target offset distance increased. This effect was more pronounced for the smaller target–applicator combinations. Selection of applicator size is based on providing coverage for the full compressed diameter of the target volume, which varies inversely with the final breast compression thickness. In all cases presented in this section, a 6 cm diameter was used as the reference size to provide coverage of the 3.5 cm target and the 8 cm applicator for the 5 cm diameter target, regardless of compression thickness.
The influence of breast compression thickness and target offset distances for a 5 cm round applicator and 3.5 cm diameter target are depicted in Figure 3. As the offset distance and compression thickness increased, the dose coverage was less uniform, but the minimum dose to the volume increased with the offset. Similarly, the dose coverage was more uniform as the offset distance and compression thickness decreased. Specifically, there were minimal differences in the 90% isodose (D90) coverage for when the target had no offset (i.e., when the target was centrally positioned as in Figure 2). Conversely, as the offset increases, the average dose to the target volume increases and D90 coverage improves.
Depictions of planar dose distributions for centralized targets are given in Figure 4. As expected, the model produced symmetric results along the beam irradiation axes. This was evident with the two-fold and four-symmetry. An additional observation in support of the model behavior was the similarity between the dose distributions presented from the CC and ML perspectives. As shown from all perspectives, the 90% isodose line included the targets for combinations of a 3.5 cm diameter target and 6 cm round applicator and a 5 cm diameter target, and 8 cm round applicator. The highest dose on the skin surface is located at the center of each beam, being 86% and 79% of the prescription dose, respectively. Also, radiation scatter increased with increasing applicator size and the dose gradients decrease.
Diameters of the D90 isodose lines for round applicators with a centrally-positioned target and uniform compressions in the CC and ML directions are shown in Table 2. With higher compression levels (i.e., thinner breast), the D90 diameters diminished in the X and Y directions. Only for the largest applicator diameter did the D90 diameters increase in the Z direction for thinner compression levels. There was X-Y symmetry for the D90 diameters with no observed differences for the largest round applicator and a maximum difference of 0.2 cm for the smallest round applicator.
The influence of non-uniform compressions between the CC and ML directions on D90 diameters is shown in Table 3 for the 6 cm round applicator with a centrally-positioned target. As in Table 2, the D90 diameters in the X and Y directions were nearly equal when the compressions in the CC and ML directions were equal. For complementary combinations of compressions in the CC and ML directions, the D90 diameters in the X and Y directions were also nearly equal. For example, with breast compression thicknesses of 3 and 8 cm in the CC and ML compressions, the D90 diameters in the X and Y directions were 4.2 and 4.8 cm, respectively, while breast compression thicknesses of 8 and 3 cm in the CC and ML compressions produced D90 diameters in the X and Y directions of 4.8 and 4.2 cm, respectively. The D90 diameters in the Z direction were largely independent of the compression asymmetries, being 4.5 ± 0.2 cm for all examined combinations of compressions.
Depictions of planar dose distributions for offset targets are given in Figure 5. While the results retained their symmetry (middle panels) from the ML perspective, dose distributions in the left and right panels exhibited asymmetry due to offsetting the target in the X direction (and subsequent repositioning of the applicators). As shown from all perspectives, the 90% isodose line included the offset targets for combinations of a 3.5 cm diameter target and 6 cm round applicator, and a 5 cm diameter target and 8 cm round applicator. For the smaller target–applicator combination, the highest dose on the skin surface was located at the center of the closest beam to the target. However, the highest dose on the skin surface for the larger target–applicator combination was located where the CC and ML beams overlap, as shown in the lower right panel. Here, the maximum skin dose was 135% of the prescription dose, though both the skin area and tissue volume associated with the overlap were small. Results from Figure 5 may be compared to Figure 4 to glean the influence of target offset, showing higher target doses with increasing offsets and smaller amounts of compression.

Effect of skin-dose optimized applicator size on target coverage and skin dose

The DVHs for covering 3.5 cm and 5 cm diameter targets with varying diameter SDO applicators are shown in Figure 6. From Figure 6A, it is evident that the smaller target is covered sufficiently by all the possible SDO applicator sizes. Figure 6B shows the underdosing that would occur when irradiating a 5 cm diameter target with a 5 cm diameter SDO applicator. Results for Figure 6 for the SDO applicator were similar to results in Figure 2 for the round applicator.
Examining target positions with offsets of 0, 2.375, and 4.75 cm from the central positions, the DVHs for covering 3.5 cm and 5.0 cm diameter targets with 6 cm and 8 cm diameter SDO applicators, respectively, are shown in Figure 6C and Figure 6D, respectively. In all cases, as with the standard round applicators, the dose coverage became less uniform as the target offset distance increased. This effect was more pronounced for the smaller target–applicator combinations. In comparison to the round applicator results in Figure 2, DVH results for the SDO applicator were more uniform.
The influence of breast compression thickness and target offset distances for a 5 cm SDO applicator and 3.5 cm diameter target are depicted in Figure 7. As the offset distance and compression thickness increases, the dose coverage is less uniform. Similarly, the dose coverage is more uniform as the offset distance and compression thickness decreased. These general behaviors follow those observed with the standard round applicators. Specifically, there are minimal differences in D90 coverage for when the target had no offset (i.e., when the target was centrally positioned as in Figure 6). Conversely, as the offset increases, the average dose to the target volume increases and D90 coverage improves. In comparison to Figure 3, DVH variations are more pronounced with SDO applicator size than for the round applicator.
Depictions of planar dose distributions for centralized targets are given in Figure 8. As shown from all perspectives, the D90 isodose line included the targets for combinations of a 3.5 cm diameter target and 6 cm SDO applicator, and a 5 cm diameter target and 8 cm SDO applicator. The highest dose on the skin surface was located adjacent to the center of each beam, being 72% and 73% of the prescription dose, respectively.
Diameters of the 90% isodose (D90) lines for SDO applicators with a centrally-positioned target and uniform compressions in the CC and ML directions are shown in Table 4. With higher compression levels (i.e., thinner breast), the D90 diameters diminished in all directions, being most prominent for the X and Y directions. These diminishments for the SDO applicators were more pronounced than for the round applicators (Table 2). There was X-Y symmetry for the D90 diameters with no observed differences for the largest round applicator and a maximum difference of 0.2 cm for the smallest SDO applicator (just like for the round applicators). These behaviors were similar to the round applicator 90% isodose line results in Table 2 for the round applicators. The D90 diameters for the SDO applicators ranged from 0.3 cm larger to 2.0 cm smaller than the D90 diameters for the round applicators, with the largest difference observed in the Z direction for the highest compression levels (i.e., 3 cm breast).
The influence of non-uniform compressions between the CC and ML directions on D90 diameters is shown in Table 5 for the 6 cm SDO applicator with a centrally-positioned target. As in Table 3, the D90 diameters in the X and Y directions were nearly equal when the compressions in the CC and ML directions were equal. For complementary combinations of compressions in the CC and ML directions, the D90 diameters in the X and Y directions were also nearly equal. For example, with breast compression thicknesses of 3 and 8 cm in the CC and ML compressions, the D90 diameters in the X and Y directions were 3.8 and 4.7 cm, respectively, while breast compression thicknesses of 8 and 3 cm in the CC and ML compressions produced D90 diameters in the X and Y directions of 4.7 and 3.6 cm, respectively. These behaviors were similar to the round applicator 90% isodose line results in Table 3 for the round applicators. Unlike for the round applicators, the D90 diameters for the SDO applicators in the Z direction were dependent of the compression asymmetries, being 3.8 ± 1.0 cm over the range of examined compression combinations. As observed in Table 4 for uniformly delivered breast compressions, the D90 diameters for the SDO applicators ranged from 0.3 cm larger to 2.0 cm smaller than the D90 diameters for the round applicators with non-uniform compressions.
Depictions of planar dose distributions for offset targets are given in Figure 9. As shown from all perspectives, the 90% isodose line included the offset targets for combinations of a 3.5 cm diameter target and 6 cm round applicator, and a 5 cm diameter target and 8 cm round applicator. For the condition having the largest offset (lower panels with X = +4.75 cm), the central aspect of the target received the lowest dose with some volumes just under D90. For the smaller target–applicator combination, the highest dose on the skin surface was located adjacent to the center of the closest beam to the target. Like for the round applicator, the highest dose on the skin surface for the larger combination of target and SDO applicator was located where the CC and ML beams overlap, as shown in the middle right panel. Here, the maximum skin dose was 130% of the prescription dose and again, the skin area and associated volumes are small. Results from Figure 9 may be compared to Figure 8 to glean the influence of target offset, generally showing higher dose hotspots on the skin and at depth near the target. Results from Figure 9 may be compared to Figure 5 to glean the influence of applicator type (i.e., internal shielding), showing the smaller irradiated volumes when using the SDO applicators.

Discussion

Interpretation of results

Using the FEA model, geometric treatment variables and their influence on composite dose distribution, target coverage, and skin dose were evaluated over a clinically meaningful range. A key finding quantified in this study was that the projected size of the target in the plane orthogonal to the treatment axis increases when the breast is under compression. This was demonstrated in the FEA model, Figure 1 and Table 1. Thus, an appropriately sized applicator must be selected to encompass the target, not in the relaxed state, but at the desired level of compression. This selection is performed by taking a low-kV mammographic radiograph with the breast in the compressed state and subsequently selecting an applicator that encompasses the target. With this approach, resultant DVHs indicate that > 90% of the target volume is generally covered by > 90% of the prescription dose, where the nominal applicator size is greater than or equal to the compressed target size. In most cases, > 95% of the target volume is covered by > 95% of the prescription dose. For APBI, an acceptable standard is > 90% of the target volume covered by > 90% of the prescription dose [48]. The only instance where dosimetric coverage was inadequate was for breast compressions  4 cm for the SDO applicator as shown in Figure 7. This phenomenon is due to the central shielding of the SDO applicators and is only manifest at extreme compressions. While the likelihood of having such thin compressions is low [27], target coverage could be made adequate by selecting the round applicator under these circumstances.
Another key finding was the dosimetric influence of target offset distance. The extent of compression had a small inverse-effect on dose coverage with the target centrally positioned, but this increased significantly with increasing lateral offset. Applicator size had only a mini­mal affect on target coverage as long as the nominal applicator size was greater than or equal to the compressed target size. Furthermore, skin dose subtended by one of the four applied fields also increased with increasing magnitude of lateral offset, though the area of the elevated dose varied substantially with the geometric conditions. Differences were minimal between the round and SDO applicators in target coverage dependence on offset distance.
Skin dose is an important consideration with the NIBB technique as a HDR 192Ir BT source (mean energy < 0.4 MeV) is used close to the skin to deliver dose at depth. Skin dose is reduced by using breast compression to limit the treatment depth and by using two orthogonal axes. To further reduce skin dose, SDO applicators were designed to augment the round applicator capabilities. These SDO applicators have a central cone design that reduces skin dose, yet maintains dose at depth by reducing skin irradiation. Evaluation of SDO applicators shows similar effect of applicator size, compression, and target lateral offset on target coverage, yet skin dose magnitude and area are reduced with SDO applicators compared to the first generation round applicators. This comes at the expense of a modest decline in SDO D90 target coverage.
When irradiating offset targets, applicator size cannot be decreased to mitigate the high skin doses resultant from applicator overlap as it is needed for target coverage. With the geometry used in the current study, an alternate beam arrangement than CC and ML delivered with all beams rotated by 45° can diminish the total surface area of irradiated skin. However, the maximum skin dose will be higher than when an applicator is positioned as close as possible to the offset target. Dose overlapped on the skin more for the round than the SDO applicators. Given these observations, clinicians should be mindful of beam orientations and target offset implications on patient skin dose.

Study limitations

A key assumption of the current model is that there is no tissue flow along the Z direction. If a more complicated strategy were implemented to add deformation along the Z direction, another assumption on such tissue properties and restricting tissue flow would be required. Given the variable tissue properties along the Z direction (such as no tissue being present beyond the breast nipple and the less-fluid characteristics as the breast couples with the chest wall), the current 2D plane-stress model was identified as the most realistic approach to the majority of the breast volume to be included in the NIBB treatment.
An additional assumption in this model is the choice for breast biomechanical properties, i.e., selected values for E and v, namely E = 30 kPa and v = 0.30. To evaluate this assumption, a separate study was conducted using alternate values of E = 45 kPa and/or v = 0.15 for the 5 cm round applicator for uniform CC and ML breast compressions of 6 cm with a 2 Gy prescription dose. Compared to the reference values of and E = 30 kPa and v = 0.30, DVH values were constant within 0.003 Gy for property values of E = 45 kPa and v = 0.30. Also, DVHs for E = 30 kPa and v = 0.15 and for E = 45 kPa and v = 0.15 were similarly equivalent, but differed from the reference values by 0.080 Gy at V95% and 0.034 Gy at V50%. Consequently, it appeared that results were generally insensitive to the selection of E value, but changed by 0.27% at V95% for choice of v value, and were generally insensitive to the selection of E value. This model behavior was expected since E dominates tissue movement in the direction of compression, while v dominates tissue flow orthogonally. The observed sensitivity of NIBB results to biomechanical properties is more favorable than what has been exhibited in other EBRT studies, where a 30% change in a biomechanical model can result in a large targeting error [48].
With EBRT, a typical challenge is the impact of shrinking tumor when performing dose accumulation across multiple treatment fractions. With NIBB, a challenge for performing dose accumulation is that the breast and target volumes may shrink upon compression, where bodily fluids exit the irradiated volume over the temporal course of the treatment. Data are needed on the constancy of the irradiated patient geometry, so that information used in advanced dose calculations may be considered robust and pertinent.
A final assumption is that the cylindrical model is representative of a true breast. The model ignores edge effects at both the chest wall, and at the distal end where the breast contour becomes more spherical rather than cylindrical. These edge effects can alter tissue motion under compression. The cylindrical geometry was employed for the current study as the true geometry and resulting edge effects vary from patient to patient. A single cylindrical diameter was assumed. Differing starting diameters will influence the resulting target size and dose distribution with compression. However, it is expected that this will minimally affect the presented results. Thus, although, breast size and shape vary from patient to patient, the employed model represents a reasonable approximation and starting point for analysis of tissue deformation and dose accumulation with the NIBB technique.

Future directions

While the current biomechanical FEA model appears to produce self-consistent results, the model requires further validation such as through direct measurements using experimental techniques, which could include patients and phantoms [49]. For example, Sioshansi et al. presented results of a clinical study tracking external fiducial markers on patients to evaluate dose summation to the skin [31]. These clinical observations of skin dose fit with benchmarking the model to the original MC results within 2%. The study of external boundaries could be extended to include tracking internal volumes such as the 3D location of in vivo fiducial clips. Validation of DIR methods has also been performed using digital or physical phantoms [50-52]. For example, gel dosimeters have been used to evaluate the accuracy of DIR methods, showing a 97% pass rate for a distance-to-agreement criterion of 0.3 cm [53].
Using incremental small steps in modeling large deformations has the drawback of introducing and accumulating numerical errors at each step. Using non-linear models that account for differences between incremental steps can address this concern, but will introduce more complexity and computational time. On the other hand, an assumption of the elasticity theory is that the initial shape of the model is stress free, which is not the case with the breast due to gravitational forces. Thus, it is desirable to remove the effect of gravity from the model to improve the simulations.
There is room for refinement of this simple model design. Instead of a uniform, tissue equivalent cylinder, CT based data can be used to create a patient-specific anatomical model that accounts for inhomogeneous material properties within the major structures within the breast. This would permit patient specific, 3D treatment planning, which could be used to refine the prescription to optimize target coverage and tissue sparing. With this future goal of patient-specific treatment planning, the interplay between model complexity and computational time would need to fit acceptable standards for clinical throughput. A first step would be to evaluate sensitivity of results to grid-size to glean the importance of high spatial resolution for the necessary standards of CT-based imaging and the required resolution for phantom validations.

Conclusions

Noninvasive image-guided breast brachytherapy is a novel method for partial breast irradiation. It has the advantage of maintaining a high degree of precision without the need for invasive catheters or applicators by using breast immobilization and high-resolution image-guidance. Furthermore, breast compression not only achieves immobilization, but also mechanically displaces non-target breast tissue out of the radiation field. To achieve conformal treatment and to reduce skin dose, treatment is delivered over two sequential nominally-orthogonal compression axes. The dose distribution in each of the compressed planes is well characterized. However, a major difficulty with this technique is the dose summation between the two treatment axes due to substantial tissue deformity under compression. Commercially-available DIR software is not able to accurately co-register imaging obtained in two states of such extreme deformation. To better understand the composite dose distribution obtained with the NIBB technique using sequential ortho­gonal compression axes for treatment, an FEA model of the breast was developed.
The FEA model developed in the current study exhibited trends in skin dose that matched MC-generated results for benchmarking within 2%, as well as observations made through clinical practice over a similar range of breast compression thicknesses and target sizes. The model provided quantitative insight on dosimetric treatment variables over a range covering of clinical circumstances. These findings highlight the need to focus on careful target localization, accurate identification of compression thickness, and target offset distance with potential dose overlap on the skin surface. As with other DIR methods, clinicians are urged to use caution if using DIR methods to determine dose accumulation for increasing clinical dose limits.

Acknowledgements

The authors express their appreciation to Ray Bricault, Jr. and Adam Brandow of Advanced Radiation Therapy, LLC for their input on this study.

Disclosure

Prof. Rivard is a shareholder of Advanced Radiation Therapy, LLC, the manufacturer of the AccuBoost® noninvasive image-guided breast brachytherapy system examin­ed herein. None of the other authors have any conflicts-of-interest for this work.

References

1. Thomadsen BR, Williamson JF, Rivard MJ et al. Anniversary paper: past and current issues, and trends in brachytherapy physics. Med Phys 2008; 35: 4708-4723.
2. Georg D, Hopfgartner J, Gòra J et al. Dosimetric considera­tions to determine the optimal technique for localized prostate cancer among external photon, proton, or carbon-ion therapy and high-dose-rate or low-dose-rate brachytherapy. Int J Radiat Oncol Biol Phys 2014; 88: 715-722.
3. Nath R, Anderson LL, Luxton G et al. Dosimetry of interstitial brachytherapy sources: recommendations of the AAPM Radiation Therapy Committee Task Group No. 43. Med Phys 1995; 22: 209-234.
4. Rivard MJ, Coursey BM, DeWerd LA et al. Update of AAPM Task Group No. 43 Report: A revised AAPM protocol for bra­chytherapy dose calculations. Med Phys 2004; 31: 633-674.
5. Rivard MJ, Beaulieu L, Mourtada F. Enhancements to commissioning techniques and quality assurance of brachytherapy treatment planning systems that use model-based dose calculation algorithms. Med Phys 2010; 37: 2645-2658.
6. Beaulieu L, Carlsson Tedgren Å, Carrier J-F et al. Report of the Task Group 186 on model-based dose calculation methods in brachytherapy beyond the TG-43 formalism: current status and recommendations for clinical implementation. Med Phys 2012; 39: 6208-6236.
7. Rivard MJ, Venselaar JLM, Beaulieu L. The evolution of bra­chytherapy treatment planning. Med Phys 2009; 36: 2136-2153.
8. Haie-Meder C, Pötter R, Van Limbergen E et al. Recommendations from Gynaecological (GYN) GEC-ESTRO Working Group (I): Concepts and terms in 3D image based 3D treatment planning in cervix cancer brachytherapy with emphasis on MRI assessment of GTV and CTV. Radiother Oncol 2005; 74: 235-245.
9. van Mourik A, van Kranen S, den Hollander S et al. Effects of setup errors and shape changes on breast radiotherapy. Int J Radiat Oncol Biol Phys 2011; 79: 1557-1564.
10. Brock KK, McShan DL, Ten Haken RK et al. Inclusion of organ deformation in dose calculations. Med Phys 2003; 30: 290-295.
11. Kirby N, Chuang C, Ueda U et al. The need for application- based adaptation of deformable image registration. Med Phys 2013; 40: 011702.
12. Kessler ML. Image registration and data fusion in radiation therapy. Brit J Radiol 2006; 79: S99-S108.
13. Pötter R, Kirisits C, Fidarova EF et al. Present status and future of high-precision image guided adaptive brachytherapy for cervix carcinoma. Acta Oncol 2008; 47: 1325-1336.
14. Tanderup K, Nesvacil N, Pötter R et al. Uncertainties in image guided adaptive cervix cancer brachytherapy: impact on planning and prescription. Radiother Oncol 2013; 107: 1-5.
15. Leibfarth S, Mönnich D, Welz S et al. A strategy for multimodal deformable image registration to integrate PET/MR into radiotherapy treatment planning. Acta Oncol 2013; 52: 1353-1359.
16. Kim H, MS, Huq MS, Houser C et al. Mapping of dose distribution from IMRT onto MRI-guided high dose rate brachytherapy using deformable image registration for cervical cancer treatments: preliminary study with commercially available software. J Contemp Brachytherapy 2014; 6: 178-184.
17. Thor M, Andersen ES, Petersen JB et al. Evaluation of an application for intensity-based deformable image registration and dose accumulation in radiotherapy. Acta Oncol 2014; 53: 1329-1336.
18. Sabater S, Andres I, Sevillano M et al. Dose accumulation during vaginal cuff brachytherapy based on rigid/deformable registration vs. single plan addition. Brachytherapy 2014; 13: 343-351.
19. Van de Kamer JB, De Leeuw AA, Moerland MA et al. Determining DVH parameters for combined external beam and brachytherapy treatment: 3D biological dose adding for patients with cervical cancer. Radiother Oncol 2010; 94: 248-253.
20. Pérez-Calatayud J, Ballester F, Das RK et al. Dose calculation for photon-emitting brachytherapy sources with average energy higher than 50 keV: Report of the AAPM and ESTRO. Med Phys 2012; 39: 2904-2929.
21. Brock KK. Results of a multi-institution deformable registration accuracy study (MIDRAS). Int J Radiat Oncol Biol Phys 2010; 76: 583-596.
22. Hepel JT, Wazer DE. A comparison of brachytherapy techniques for partial breast irradiation. Brachytherapy 2012; 11: 163-175.
23. Hepel JT, Hiatt JR, Sha S et al. The rationale, technique, and feasibility of partial breast irradiation using noninvasive image-guided breast brachytherapy. Brachytherapy 2014; 13: 493-501.
24. Yang Y, Rivard MJ. Monte Carlo simulations and radiation dosimetry measurements of peripherally applied HDR 192Ir breast brachytherapy D-shaped applicators. Med Phys 2009; 36: 809-815.
25. Rivard MJ, Melhus CS, Wazer DE, Bricault Jr RJ. Dosimetric characterization of round HDR 192Ir AccuBoost applicators for breast brachytherapy. Med Phys 2009; 36: 5027-5032.
26. Yang Y, Rivard MJ. Dosimetric optimization of a conical breast brachytherapy applicator for improved skin-dose sparing. Med Phys 2010; 37: 5665-5671.
27. Hamid S, Rocchio K, Arthur D et al. A multi-institutional study of feasibility, implementation, and early clinical results with noninvasive breast brachytherapy for tumor bed boost. Int J Radiat Oncol Biol Phys 2012; 83: 1374-1380.
28. Leonard KL, Hepel JT, Styczynski JR et al. Breast boost using noninvasive image-guided breast brachytherapy vs. external beam: a 2:1 matched-pair analysis. Clin Breast Cancer 2013; 13: 455-459.
29. Yang Y, Melhus CS, Sioshansi S et al. Treatment planning of a skin-sparing conical breast brachytherapy applicator using conventional brachytherapy software. Med Phys 2011; 38: 1519-1525.
30. Rivard MJ, Melhus CS, Granero D et al. An approach to using conventional brachytherapy software for clinical treatment planning of complex, Monte Carlo-based brachytherapy dose distributions. Med Phys 2009; 36: 1968-1975.
31. Sioshansi S, Rivard MJ, Hiatt JR et al. Dose modeling of noninvasive image-guided breast brachytherapy in comparison to electron beam boost and three-dimensional conformal accelerated partial breast irradiation. Int J Radiat Oncol Biol Phys 2011; 80: 410-416.
32. Schnabel J, Tanner C, Castellano-Smith A, et al. Validation of nonrigid image registration using finite-element methods: application to breast MR images. IEEE Trans Med Imaging 2003; 22: 238-247.
33. Chung JH, Rajagopal V, Nielsen PMF, Nash M. A biomechanical model of mammographic compression. Biomech Model Mechanobiol 2007; 7: 43-52.
34. del Palomar P, Calvo A, Herrero B et al. A finite element model to accurately predict real deformations of the breast. Med Eng Phys 2008; 30: 1089-1097.
35. Azar F, Metaxas D, Schnall M. A finite element model of the breast for predicting mechanical deformations during biopsy procedures. In: Mathematical Methods in Biomedical Image Analysis Proceedings, 2000, pp. 38-45.
36. Samani A, Bishop J, Yaffe M, Plewes D. Biomedical 3D finite element modeling of the human breast using MRI data. IEEE Trans Med Imaging 2001; 20: 271-279.
37. Rajagopal V, Lee A. Creating individual-specific biomechanical models of the breast for medical image analysis. Acad Radiol 2008; 15: 1425-1436.
38. Krouskop TA, Wheeler TM, Kallel F et al. Elastic moduli of breast and prostate tissues under compression. Ultrasonic Imaging 1998; 20: 260-274.
39. Bakic PR. Breast tissue description and modeling in mammography. PhD thesis, Lehigh University, 2000.
40. Kuhlmanna M, Fear EC. Mechanical model of the breast for the prediction of deformation during imaging. Med Engin Phys 2013; 35: 470-478.
41. Hu T, Desai JP. Soft-tissue material properties under large deformation: strain rate effect. Conf Proc IEEE Eng Med Biol Soc 2004; 4: 2758-2761.
42. Li M, Castillo E, Zheng XL et al. Modeling lung deformation: a combined deformable image registration method with spatially varying Young’s modulus estimates. Med Phys 2013; 40: 081902.
43. Shewchuk JR. Delaunay refinement algorithms for triangular mesh generation. Comp Geom: Theory Appl 2002; 22: 21-74.
44. Plewes D, Bishop J. Visualization and quantification of breast cancer biomechanical properties with magnetic resonance elastography. Phys Med Biol 2000; 45: 1591-1610.
45. Devi C, Chandran R. Measurement of visco-elastic properties of breast-tissue mimicking materials using diffusing wave spectroscopy. J Biomed Optics 2007; 12: 034035.
46. Van Houten E, Doyley M, Kennedy FE et al. Initial in vivo experience with steady-state subzone-based MR elastography of the human breast. J Magn Res Imag 2003; 17: 72-85.
47. Fung YC. A First Course in Continuum Mechanics. Prentice Hall, Englewood Cliffs 1994.
48. Chi Y, Liang J, Yan D. A material sensitivity study on the accuracy of deformable organ registration using linear biomechanical models. Med Phys 2006; 33: 421-433.
49. Hoffmann C, Krause S, Stoiber EM et al. Accuracy quantification of a deformable image registration tool applied in a clinical setting. J Appl Clin Med Phys 2014; 15: 4564.
50. Kashani R, Hub M, Balter JM et al. Objective assessment of deformable image registration in radiotherapy: A multi-institution study. Med Phys 2008; 35: 5944-5953.
51. Serban M, Heath E, Stroian G et al. A deformable phantom for 4D radiotherapy verification: design and image registration evaluation. Med Phys 2008; 35: 1094-1102.
52. Hoffmann C, Krause S, Stoiber EM et al. Accuracy quantification of a deformable image registration tool applied in a clinical setting. J Appl Clin Med Phys 2014; 15: 4564.
53. Vicini F, Beitsch PD, Quiet CA et al. Three-year analysis of treatment efficacy, cosmesis, and toxicity by the American Society of Breast Surgeons MammoSite Breast Brachytherapy Registry Trial in patients treated with accelerated partial breast irradiation (APBI). Cancer 2008; 112: 758-766.