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

PURPOSE: Noninvasive image-guided breast brachytherapy delivers conformal HDR (192)Ir brachytherapy treatments with the breast compressed, and treated in the cranial-caudal and medial-lateral directions. This technique subjects breast tissue to extreme deformations not observed for other disease sites. Given that, commercially-available software for deformable image registration cannot accurately co-register image sets obtained in these two states, a finite element analysis based on a biomechanical model was developed to deform dose distributions for each compression circumstance for dose summation.
MATERIAL AND METHODS: The model assumed the breast was under planar stress with values of 30 kPa for Young's modulus and 0.3 for Poisson's ratio. Dose distributions from round and skin-dose optimized applicators in cranial-caudal and medial-lateral compressions were deformed using 0.1 cm planar resolution. Dose distributions, skin doses, and dose-volume histograms were generated. Results were examined as a function of breast thickness, applicator size, target size, and offset distance from the center.
RESULTS: Over the range of examined thicknesses, target size increased several millimeters as compression thickness decreased. This trend increased with increasing offset distances. Applicator size minimally affected target coverage, until applicator size was less than the compressed target size. In all cases, with an applicator larger or equal to the compressed target size, > 90% of the target covered by > 90% of the prescription dose. In all cases, dose coverage became less uniform as offset distance increased and average dose increased. This effect was more pronounced for smaller target-applicator combinations.
CONCLUSIONS: The model exhibited skin dose trends that matched MC-generated benchmarking results within 2% and clinical observations over a similar range of breast thicknesses and target sizes. The model provided quantitative insight on dosimetric treatment variables over a range of clinical circumstances. These findings highlight the need for careful target localization and accurate identification of compression thickness and target offset.

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 distributions 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 Symphony™ by MIM Software, Inc. (Cleveland, OH, USA), Velocity™ by Varian Medical Systems, Inc. (Palo Alto, CA, USA), RTx™ by Mirada Medical (Oxford, UK), and IMSimQA™ 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 ranked 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:

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 MATLAB™ (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 [] = [K] · [], where is nodal force vector and 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.

Fig. 1

A) Calculation geometry illustrating the coordinate system and reference axes. The front circular surface depicts a mesh of 1,300 triangular elements representing a 15 cm diameter circular model of 2D plane-stress analysis with 0.1 cm planar thickness. B) Parallel-opposed beams from round AccuBoost applicators are applied along the Y-axis when the cylindrical breast model is compressed. C) The breast model is relaxed to the cylindrical shape then compressed in the X-axis direction. D) Parallel-opposed beams from round AccuBoost applicators are applied along the X-axis. E) The breast model is relaxed to the reference cylindrical shape, depicting dose summation from the applicators delivered along dual orthogonal vectors. It is evident that the projected target size increases with compression when comparing part E) with parts B) and C)

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 . 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 orthogonal 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.

Table 1

Observed diameters in the cranial-caudal (CC) and medial-lateral (ML) views of compressed targets given varying initial (uncompressed) diameters, compressed breast thicknesses, and offset distances. All units are in centimeters

Target diameter	Breast compression	Offset = 0 cm (centered target)		Offset = 2.375 cm (X-direction)		Offset = 4.75 cm (X-direction)
		CC	ML	CC	ML	CC	ML
2	3	2.80	2.82	2.87	2.69	2.77	2.32
	4	2.79	2.81	2.84	2.68	2.57	2.32
	5	2.77	2.79	2.78	2.67	2.41	2.31
	6	2.73	2.75	2.70	2.63	2.29	2.29
	7	2.67	2.68	2.60	2.58	2.20	2.27
	8	2.59	2.60	2.49	2.52	2.14	2.25
3.5	3	4.98	4.93	5.08	4.80	5.01	4.35
	4	4.96	4.91	5.02	4.79	4.77	4.34
	5	4.92	4.87	4.93	4.75	4.53	4.32
	6	4.84	4.79	4.79	4.68	4.31	4.28
	7	4.73	4.68	4.62	4.58	4.13	4.22
	8	4.58	4.53	4.44	4.45	3.99	4.14
5	3	7.12	7.13	7.25	6.88	7.17	6.46
	4	7.09	7.10	7.14	6.85	6.91	6.44
	5	7.02	7.02	6.97	6.79	6.63	6.40
	6	6.90	6.90	6.75	6.69	6.37	6.33
	7	6.72	6.72	6.51	6.54	6.12	6.22
	8	6.50	6.50	6.26	6.35	5.89	6.09

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 cm3) 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.

Fig. 2

Comparative dose value histograms (DVHs) showing the influence of varying the diameter of round applicators on coverage for centrally-positioned spherical targets having A) 3.5 cm and B) 5 cm diameters when uncompressed, and also for targets offset by 0, 2.375, and 4.75 cm from the central position for C) 3.5 cm diameter target with a 6 cm diameter round applicator and D) 5 cm diameter target with a 8 cm diameter round applicator. In all cases, the breast was compressed to a thickness of 6 cm

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 (D ) 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 D coverage improves.

Fig. 3

Comparative dose value histograms (DVHs) showing the influence of breast thickness and target offset distances of A) 0 cm, B) 2.375 cm, and C) 4.75 cm for a 5 cm round applicator and 3.5 cm diameter target

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.

Fig. 4

For centralized targets (dark blue circles) irradiated with four beams and with 6 cm of compression, the 2D dose distributions are depicted (upper panels) for a 3.5 cm diameter target and the 6 cm round applicator, and (lower panels) for a 5 cm diameter target and the 8 cm round applicator. Dose distributions on the surface are depicted on the right panels. Results on the XZ plane were not significantly different from results on the YZ plane. Theta = 0° is in the positive X direction

Diameters of the D 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 D diameters diminished in the X and Y directions. Only for the largest applicator diameter did the D diameters increase in the Z direction for thinner compression levels. There was X-Y symmetry for the D diameters with no observed differences for the largest round applicator and a maximum difference of 0.2 cm for the smallest round applicator.

Table 2

Diameters of the 90% isodose lines for round applicators with a centrally-positioned target and uniform compressions in the cranial-caudal and medial-lateral directions. All units are in centimeters

Breast compression	5 round			6 round			7 round			8 round
	X	Y	Z	X	Y	Z	X	Y	Z	X	Y	Z
3	3.6	3.6	3.9	4.1	4.0	4.7	4.9	4.9	5.5	5.6	5.6	6.7
4	3.6	3.6	3.9	4.1	4.0	4.3	5.0	4.9	5.1	5.6	5.6	6.3
5	3.6	3.8	3.9	4.3	4.2	4.5	5.0	5.0	5.1	5.8	5.8	6.3
6	4.0	4.0	4.3	4.6	4.6	4.5	5.2	5.2	5.3	6.2	6.2	6.1
7	4.2	4.1	4.1	4.9	4.8	4.7	5.6	5.6	4.9	6.4	6.4	5.7
8	4.6	4.4	4.3	5.2	5.2	4.5	6.1	6.1	5.3	7.0	7.0	6.1

The influence of non-uniform compressions between the CC and ML directions on D diameters is shown in Table 3 for the 6 cm round applicator with a centrally positioned target. As in Table 2, the D 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 D 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 D 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 D diameters in the X and Y directions of 4.8 and 4.2 cm, respectively. The D diameters in the Z direction were largely independent of the compression asymmetries, being 4.5 ± 0.2 cm for all examined combinations of compressions.

Table 3

Diameters of the 90% isodose lines for the 6 cm round applicator with a centrally-positioned target and non-uniform compressions in the cranial-caudal (CC) and medial-lateral (ML) directions. All units are in centimeters

Breast compression	3 ML			4 ML			5 ML			6 ML			7 ML			8 ML
	X	Y	Z	X	Y	Z	X	Y	Z	X	Y	Z	X	Y	Z	X	Y	Z
3 CC	4.1	4.0	4.7	4.1	4.0	4.5	4.2	4.2	4.5	4.2	4.6	4.5	4.2	4.6	4.7	4.2	4.8	4.5
4 CC	4.2	4.0	4.5	4.1	4.0	4.3	4.2	4.2	4.3	4.2	4.6	4.5	4.3	4.6	4.5	4.2	4.8	4.5
5 CC	4.2	4.0	4.5	4.2	4.0	4.3	4.3	4.2	4.5	4.4	4.6	4.5	4.5	4.7	4.5	4.5	4.9	4.5
6 CC	4.6	4.1	4.5	4.5	4.1	4.5	4.6	4.4	4.5	4.6	4.6	4.5	4.8	4.8	4.5	4.8	5.0	4.5
7 CC	4.6	4.2	4.7	4.6	4.2	4.5	4.8	4.5	4.5	4.8	4.7	4.5	4.9	4.8	4.7	4.9	5.1	4.5
8 CC	4.8	4.2	4.5	4.8	4.2	4.5	4.9	4.5	4.5	5.0	4.8	4.5	5.2	4.9	4.5	5.2	5.2	4.5

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.

Fig. 5

For X-offset targets (dark blue circles) irradiated with four beams and with 6 cm of compression, the 2D dose distributions are depicted (upper panels) for a 3.5 cm diameter target (+2.38 cm offset) and the 6 cm round applicator, and (lower panels) for a 5 cm diameter target (+2.51 cm offset) and the 8 cm round applicator. Dose distributions on the surface are depicted on the right panels. Results on the XZ plane were not significantly different from results on the YZ plane. Theta = 0° is in the positive X direction

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.

Fig. 6

Comparative dose value histograms (DVHs) showing the influence of varying the diameter of skin-dose optimized (SDO) applicators on coverage for centrally-positioned spherical targets having A) 3.5 cm and B) 5 cm diameters when uncompressed, and also for targets offset by 0, 2.375, and 4.75 cm from the central position for C) 3.5 cm diameter target with a 6 cm diameter SDO applicator and D) 5 cm diameter target with a 8 cm diameter SDO applicator. In all cases, the breast was compressed to a thickness of 6 cm

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 D 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 D coverage improves. In comparison to Figure 3, DVH variations are more pronounced with SDO applicator size than for the round applicator.

Fig. 7

Comparative dose value histograms (DVHs) showing the influence of breast thickness and target offset distances of A) 0 cm, B) 2.375 cm, and C) 4.75 cm for a 5 cm SDO applicator and 3.5 cm diameter target

Depictions of planar dose distributions for centralized targets are given in Figure 8. As shown from all perspectives, the D 90 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.

Fig. 8

For centralized targets (dark blue circles) irradiated with four beams and with 6 cm of compression, the 2D dose distributions are depicted (upper panels) for a 3.5 cm diameter target and the 6 cm skin-dose optimized (SDO) applicator, and (lower panels) for a 5 cm diameter target and the 8 cm SDO applicator. Dose distributions on the surface are depicted on the right panels. Results on the XZ plane were not significantly different from results on the YZ plane. Theta = 0° is in the positive X direction

Diameters of the 90% isodose (D ) 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 D 90 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 D 90 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 D diameters for the SDO applicators ranged from 0.3 cm larger to 2.0 cm smaller than the D diameters for the round applicators, with the largest difference observed in the Z direction for the highest compression levels (i.e., 3 cm breast).

Table 4

Diameters of the 90% isodose lines for skin-dose optimized (SDO) applicators with a centrally-positioned target and uniform compressions in the cranial-caudal and medial-lateral directions. All units are in centimeters

Breast compression	5 SDO			6 SDO			7 SDO			8 SDO
	X	Y	Z	X	Y	Z	X	Y	Z	X	Y	Z
3	2.4	2.4	2.1	3.0	3.0	2.7	4.2	4.2	3.5	5.3	5.2	4.7
4	3.0	3.0	3.1	3.6	3.6	3.7	4.4	4.4	4.5	5.4	5.4	5.3
5	3.4	3.4	3.5	3.7	3.8	3.9	4.8	4.9	5.1	5.5	5.5	5.9
6	3.6	3.7	3.7	4.4	4.4	4.3	5.0	4.9	4.7	6.0	5.9	5.9
7	4.0	4.0	3.9	4.8	4.9	4.5	5.6	5.6	4.9	6.4	6.2	5.7
8	4.6	4.6	4.1	5.4	5.5	4.7	6.0	6.0	4.9	7.0	7.0	5.7

The influence of non-uniform compressions between the CC and ML directions on D diameters is shown in Table 5 for the 6 cm SDO applicator with a centrally-positioned target. As in Table 3, the D 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 D 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 D 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 D 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 D 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 D diameters for the SDO applicators ranged from 0.3 cm larger to 2.0 cm smaller than the D diameters for the round applicators with non-uniform compressions.

Table 5

Diameters of the 90% isodose lines for the 6 cm skin-dose optimized (SDO) applicator with a centrallypositioned target and non-uniform compressions in the cranial-caudal (CC) and medial-lateral (ML) directions. All units are in centimeters

Breast compression	3 ML			4 ML			5 ML			6 ML			7 ML			8 ML
	X	Y	Z	X	Y	Z	X	Y	Z	X	Y	Z	X	Y	Z	X	Y	Z
3 CC	3.0	3.0	2.7	3.2	3.4	3.1	3.0	3.6	3.1	3.4	4.0	3.1	3.6	4.5	3.3	3.8	4.7	3.3
4 CC	3.4	3.4	3.1	3.6	3.6	3.7	3.4	3.8	3.7	3.6	4.2	3.9	3.8	4.6	3.9	4.0	4.9	3.9
5 CC	3.6	3.0	3.1	3.8	3.4	3.7	3.7	3.8	3.9	3.9	4.2	4.1	4.1	4.6	4.1	4.2	4.9	4.3
6 CC	4.1	3.4	3.1	4.3	3.6	3.9	4.2	3.8	4.1	4.4	4.4	4.3	4.6	4.8	4.3	4.6	5.1	4.5
7 CC	4.6	3.6	3.3	4.6	3.8	3.9	4.6	4.0	4.1	4.8	4.6	4.3	4.8	4.9	4.5	5.0	5.3	4.5
8 CC	4.7	3.6	3.3	5.0	3.8	3.9	5.0	4.2	4.3	5.2	4.6	4.5	5.3	5.0	4.5	5.4	5.5	4.7

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 D . 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.

Fig. 9

For X-offset targets (dark blue circles) irradiated with four beams and with 6 cm of compression, the 2D dose distributions are depicted (upper panels) for a 3.5 cm diameter target (+2.38 cm offset) and the 6 cm skin-dose optimized (SDO) applicator, (middle panels) for a 5 cm diameter target (+2.51 cm offset) and the 8 cm SDO applicator, and (lower panels) for a 3.5 cm diameter target (+4.75 cm offset) and the 6 cm SDO applicator. Dose distributions on the surface are depicted on the right panels. Dose distributions on the surface are depicted on the right panels. Results on the XZ plane were not significantly different from results on the YZ plane. Theta = 0° is in the positive X direction

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 minimal 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 D 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 V and 0.034 Gy at V . Consequently, it appeared that results were generally insensitive to the selection of E value, but changed by 0.27% at V 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 orthogonal 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.