Brachy 2D Algorithm
The Brachy 2D web application is designed to teach basic concepts about HDR brachytherapy planning by providing a simple, adjustable test case. To help further understand how brachytherapy treatment planning works, this article describes how the app takes user input and calculates the corresponding dose distribution. The first section describes the system geometry and input variables, followed by the second section which describes the dose calculation itself in more detail.
I. System Geometry & User Input
The system geometry is defined using a two-dimensional coordinate system. The x-axis and y-axis are horizontal and vertical, respectively, with the origin (0, 0) located at the center of the image. All distance units are given as centimeters (cm). As explained on the web app page itself, the interstitial needle has 9 dwell positions. These positions are located along the y-axis (y = 0) with x ranging from -5 to -1, at 0.5 cm intervals. In addition to the target/OAR design, the user also inputs the prescription dose D_{Rx} and the prescription point (-3, d_{Rx}). Note that -3 is constant while d_{Rx} can be adjusted. Finally, the relative dwell weights w_i are also assigned by the user.
II. Dose Calculation
The Brachy 2D web app calculates the dose from brachytherapy sources using the AAPM TG-43 formalism. Using this formalism, the dose rate \dot{D} from a source can be calculated at any point defined by a radius r and an angle \theta using the following formula:
\dot{D}(r,\theta) = S_K \Lambda \frac{G(r,\theta)}{G(r_0,\theta_0)} g(r) F(r,\theta)
where:
r: radial distance from the calculation point to the source center
\theta: angle between r and the source's central axis
S_K: Strength of the brachytherapy source, in terms of air kerma strength
\Lambda: Dose rate constant
\frac{G(r,\theta)}{G(r_0,\theta_0)}: Geometry factor
g(r): Radial dose function
F(r,\theta): Anisotropy function
Please refer to the AAPM TG-43 reports for further explanation of these terms. The air kerma strength is set as a constant in the web app, with a value of 15,000 cGy cm2 hr-1. The TG-43 equation is implemented by the BrachyDose43 code, which automatically calculates the dose rate for a given source provided S_K, r, and \theta. Please refer to the source code for more information on BrachyDose43 calculation. Therefore the remaining task for the app is to provide values of r and \theta, and then the dose rate from a single dwell position at that point can be calculated. While polar coordinates are used in the TG-43 equation, the calculation and prescription points for the app are in Cartesian coordinates. Assuming that the source travels along the x-axis, with dwell positions defined as (s_i, 0) for 1 \le i \le 9, r and \theta can be calculated for any Cartesian point (x, y) as follows:
r_i(x, y) = \sqrt{(x-s_i)^2 + y^2}
\theta_i(x, y) = cos^{-1} \big\{ (x-s_i) / r_i \big\}
A. Dose Normalization
The TG-43 equation calculates the dose rate to a point, but we actually want the dose. Therefore, before calculating the 2D dose distribution the actual dwell times, with the appropriate weights, need to be calculated. This is the purpose of the prescription dose and point; the dwell times will be calculated so that the dose at the prescription point is equal to the prescription dose. Since the dwell weights w_i are relative, they will first need to be scaled by their sum W, which is calculated as W = \sum^{9}_{i=1} w_i. The dose rate at the prescription point for the given dwell weight configuration can be calculated as the sum of contributions from all dwell positions, using the TG-43 equation:
\dot{D}_{Rx} = \sum^9_{i=1} \dot{D}(r_i[-3, d_{Rx}], \theta_i[-3, d_{Rx}])
Given this dose rate we can calculate the total time T needed to deliver the prescription dose: T = D_{Rx} / \dot{D}_{Rx}. This time is distributed according to the dwell weights, so that the dwell time for position i, t_i, is equal to (w_iT)/W
B. 2D Dose Calculation
Having obtained the actual dwell times based on the prescription point, the 2D dose can be calculated. First, a 2D dose grid needs to be defined. The dose grid for this app has dimensions of 240 \times 240 pixels, with each pixel having dimensions of 1 mm \times 1 mm. This dose grid is centered at the origin. For the m_{th} row and n_{th} column, the Cartesian coordinates x and y of a pixel's center can be calculated accordingly:
x_n = -N_x / 2 + p_x / 2 + np_x
y_m = -N_y / 2 + p_y / 2 + mp_y
Using the same assumptions as above, r and \theta can be calculated for any dose grid pixel as follows:
r_{mn}(i) = \sqrt{[x_n-d_x(i)]^2 + [y_m-d_y(i)]^2}
\theta_{mn}(i) = cos^{-1} \big( [x_n-d_x(i)] / r_{mn}(i) \big)
The values r_{mn}(i) and \theta_{mn}(i) can be used directly in the TG-43 dose calculation equation, and adding the dose rates at a point from all dwell positions results in the total dose rate:
D(r_{mn}, \theta_{mn}) = \sum^{9}_{i=1} t_i \dot{D}_i(r_{mn}(i), \theta_{mn}(i))
This process is repeated for each point to produce the 2D dose grid.