# Finite-element modeling of the alpha particle dose of realistic sources used
in Diffusing Alpha-emitters Radiation Therapy

G. Heger, L. Arazi111Corresponding Author. Unit of Nuclear Engineering,
Faculty of Engineering Sciences, Ben-Gurion University of the Negav, P.O.B.
653 Be’er-Sheva 8410501, Israel<EMAIL_ADDRESS>
###### Abstract

Diffusing Alpha-emitters Radiation Therapy (DaRT) is a new method for treating
solid tumors using alpha particles. Unlike conventional radiotherapy, where
the physical models for dose calculations are known and routinely used, for
DaRT a new framework, called the Diffusion-Leakage (DL) model, had to be
developed. In this work we provide a detailed description of a finite-element
numerical scheme for solving the time-dependent DL model equations for
cylindrical DaRT sources (“seeds”) of finite diameter and length in two
dimensions. Using a fully-implicit scheme and adaptive time step, this
approach allows to accurately follow temporal transients ranging from seconds
to days. In addition to the full two-dimensional calculation, we further
provide a closed-form approximation and a simple one-dimensional scheme to
solve the equations for infinitely-long cylindrical sources. These simpler
solutions can be used both to validate the two-dimensional code and to enable
efficient parameter scans to study the properties of DaRT seed lattices. We
compare these approximations to the full 2D solution over the relevant
parameter space, providing guidelines on their usability and limitations.

††: Phys. Med. Biol.

  * October 2021

##  1 Introduction

Diffusing Alpha-Emitters Radiation Therapy (DaRT), first introduced by ?, is a
new form of radiation therapy, which enables the treatment of solid tumors
using alpha particles. In DaRT, tumors are treated with sources (“seeds”)
carrying low activities of 224Ra below their surface. Once inside the tumor,
the seeds continuously release from their surface the short-lived daughter
atoms of 224Ra: 220Rn, 216Po, 212Pb, 212Bi, 212Po and 208Tl. These spread by
diffusion (with possible contribution by convective effects), creating –
primarily through their alpha decays – a “kill-region” measuring several mm in
diameter around each seed. The full decay chain, along with the isotopes’
half-lives, decay modes and mean alpha particle energies is shown in figure 1.

Figure 1: The 224Ra decay chain. Data taken from the NuDat2 database website,
https://www.nndc.bnl.gov/nudat2/

DaRT was, and still is, extensively investigated in in vitro and in vivo
preclinical studies on a large number of cancer types, as a stand-alone
treatment (?, ?, ?, ?, ?), in combination with chemotherapy (?, ?, ?, ?), and
as a stimulator of local and systemic anti-tumor immune response (?, ?, ?, ?,
?, ?, ?). Since 2017, DaRT is under clinical investigations in human patients,
starting with locally advanced and recurrent squamous cell carcinoma of the
skin and head and neck (?). Results of the first-in-man trial were highly
promising in terms of both efficacy and safety: all treated tumors shrank
drastically (by $\sim 30-100\%$), beginning in the first few days after the
treatment, with $\sim 80\%$ of the tumors exhibiting complete response.
Adverse affects were mild to moderate, with no observable local or systemic
radiation-induced damage to healthy tissue. The alpha dose to all organs,
resulting from 212Pb leaving the tumor through the blood, was calculated to be
on the cGy level, with blood and urine activity measurements consistent with
the prediction of a biokinetic model (?). In one patient, untreated lesions
shrank and disappeared when another was treated with DaRT (in the absence of
any other treatment), suggesting an abscopal effect (?).

An in-depth description of the underlying physics of DaRT was given in (?). As
a first step towards treatment planning, it introduced a simplified
theoretical approach: the “Diffusion-Leakage (DL) model”. Its underlying
assumptions are that the medium into which the daughters of 224Ra are released
is homogeneous, isotropic and does not change with time, and that convective
effects have a short correlation length and can therefore be described by an
effective isotropic diffusion term. A further assumption is that 212Pb and
212Bi can be cleared from the tumor through the blood at a uniform local rate.
The discussion in (?) showed that of the six radionuclides released from the
source, one only needs to model the migration of 220Rn, 212Pb and 212Bi as
their respective short-lived daughters are in local secular equilibrium.

The focus of (?) was the alpha particle dose of the ideal point source, which
constitutes the basic building block for any arbitrary configuration of line
sources. In this work, we expand the treatment to realistic cylindrical seeds
of finite diameter and length. We derive closed-form asymptotic and
approximate time-dependent solutions to the diffusion-leakage model for
infinite cylindrical sources, and describe in detail numerical schemes for
solving the complete time-dependent problem in both one and two dimensions, to
model the alpha particle dose for realistic seed geometries. We compare the 2D
numerical solution to the 1D case and to other approximations, and provide
examples for calculating the alpha particle dose in hexagonal seed lattices.
The presented calculations can serve for generating dose lookup tables as a
pragmatic approach to treatment planning.

##  2 The diffusion-leakage model in cylindrical coordinates

The underlying assumptions of the DL model, described at length in (?), are:

  * •

The migration of atoms inside the tumor is governed by diffusion.

  * •

The tissue is homogeneous, isotropic, and time-independent, and thus the
diffusion coefficients are constant.

  * •

It is sufficient to model the migration of 220Rn, 212Pb and 212Bi. 216Po,
212Po and 208Tl are in local secular equilibrium with their respective parent
isotopes (with suitable branching ratios for the latter two).

  * •

212Pb migration can be described using a single effective diffusion
coefficient.

  * •

212Pb atoms reaching major blood vessels are trapped in red blood cells and
are then immediately cleared from the tumor. This is represented by a single
constant sink term.

  * •

220Rn atoms do not form chemical bonds and are very short-lived, and therefore
the equation for Rn does not include a sink term.

  * •

The equation for 212Bi does contain a sink term, but it is considered a
second-order effect and generally set to zero.

We consider the case of a cylindrical source of radius $R_{0}$ and length $l$
along the $z$ axis, and assume axial symmetry. Under the above assumptions, in
cylindrical coordinates $(r,z)$ the equations describing the dynamics of the
main daughter atoms in the decay chain – 220Rn, 212Pb and 212Bi – are:

$\qquad\qquad\frac{\partial n_{Rn}}{\partial
t}=D_{Rn}\left(\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial
n_{Rn}}{\partial r}\right)+\frac{\partial^{2}n_{Rn}}{\partial
z^{2}}\right)-\lambda_{Rn}n_{Rn}$ (1)

$\qquad\qquad\frac{\partial n_{Pb}}{\partial
t}=D_{Pb}\left(\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial
n_{Pb}}{\partial r}\right)+\frac{\partial^{2}n_{Pb}}{\partial
z^{2}}\right)+\lambda_{Rn}n_{Rn}-\left(\lambda_{Pb}+\alpha_{Pb}\right)n_{Pb}$
(2)

$\qquad\qquad\frac{\partial n_{Bi}}{\partial
t}=D_{Bi}\left(\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial
n_{Bi}}{\partial r}\right)+\frac{\partial^{2}n_{Bi}}{\partial
z^{2}}\right)+\lambda_{Pb}n_{Pb}-\left(\lambda_{Bi}+\alpha_{Bi}\right)n_{Bi}$
(3)

where $n_{Rn}$, $n_{Pb}$, $n_{Bi}$, $D_{Rn}$, $D_{Pb}$, $D_{Bi}$ and
$\lambda_{Rn}$, $\lambda_{Pb}$, $\lambda_{Bi}$ are the number densities,
diffusion coefficients and decay rate constants of 220Rn, 212Pb and 212Bi,
respectively. $\alpha_{Pb}$ and $\alpha_{Bi}$ are the leakage rate constants,
accounting for clearance through the blood, of 212Pb and 212Bi. As discussed
in (?), one can generally assume $\alpha_{Bi}=0$, as will be done here. The
boundary conditions, for $r\rightarrow R_{0}$ and $|z|\leq l/2$ ($z=0$ at the
seed mid plane), are:

$\lim_{r\rightarrow R_{0}}2\pi
rj_{Rn}\left(r,z,t\right)=P_{des}(Rn)\frac{\Gamma_{Ra}^{src}(0)}{l}e^{-\lambda_{Ra}t}$
(4)

$\lim_{r\rightarrow R_{0}}2\pi
rj_{Pb}(r,z,t)=\left(P_{des}^{eff}(Pb)-P_{des}(Rn)\right)\frac{\Gamma_{Ra}^{src}(0)}{l}e^{-\lambda_{Ra}t}$
(5)

$\lim_{r\rightarrow R_{0}}2\pi rj_{Bi}(r,z,t)=0$ (6)

where $j_{x}=-D_{x}\partial{n_{x}}/\partial{r}$, with $x$ representing 220Rn,
212Pb and 212Bi. $\Gamma_{Ra}^{src}(0)$ is the initial 224Ra activity on the
source (assumed to be uniform) and $\lambda_{Ra}$ is the 224Ra decay rate
constant. $P_{des}(Rn)$ and $P_{des}^{eff}(Pb)$ are the desorption
probabilities of 220Rn and 212Pb, respectively, representing the probability
that a decay of a 224Ra on the source will lead to the emission of either a
220Rn or 212Pb atom into the tumor; we use the term “effective” for
$P_{des}^{eff}(Pb)$, because it includes several emission pathways, as
discussed in (?). For $|z|>l/2$, $\lim_{r\rightarrow R_{0}}rj_{x}(r,z,t)=0$
for the three isotopes.

The solution for equations (1-3) provides the number densities
$n_{Rn}(r,z,t)$, $n_{Pb}(r,z,t)$ and $n_{Bi}(r,z,t)$. The alpha dose is
calculated under the assumption that the range of alpha particles is much
smaller than the dominant diffusion length of the problem (see below). The
dose from source insertion to time $t$ is comprised of two contributions: one
arising from the summed alpha particle energy of the pair 220Rn+216Po, and the
other from the the alpha decay of either 212Bi or 212Po:

$Dose_{\alpha}(RnPo;\,r,z,t)=\frac{E_{\alpha}(RnPo)}{\rho}\int_{0}^{t}\lambda_{Rn}n_{Rn}(r,z,t^{\prime})dt^{\prime}$
(7)

$Dose_{\alpha}(BiPo;\,r,z,t)=\frac{E_{\alpha}(BiPo)}{\rho}\int_{0}^{t}\lambda_{Bi}n_{Bi}(r,z,t^{\prime})dt^{\prime}$
(8)

where $E_{\alpha}(RnPo)=(6.288+6.778)\>\textrm{MeV}=13.066\>\textrm{MeV}$ is
the total alpha particle energy of 220Rn and 216Po,
$E_{\alpha}(BiPo)=7.804\>\textrm{MeV}$ is the weighted-average energy of the
alpha particles emitted by 212Bi and 212Po, and $\rho$ is the tissue density.
In what follows, we define the “asymptotic dose” as the dose delivered from
source insertion to infinity – in practice, over several weeks.

The discussion in (?) has shown the the spread of 220Rn, 212Pb and 212Bi is
governed by their respective diffusion lengths, defined as:

$L_{Rn}=\sqrt{\frac{D_{Rn}}{\lambda_{Rn}-\lambda_{Ra}}}$ (9)

$L_{Pb}=\sqrt{\frac{D_{Pb}}{\lambda_{Pb}+\alpha_{Pb}-\lambda_{Ra}}}$ (10)

$L_{Bi}=\sqrt{\frac{D_{Bi}}{\lambda_{Bi}+\alpha_{Bi}-\lambda_{Ra}}}$ (11)

The diffusion length is a gross measure of the average displacement of an atom
from its creation site, to the point of its decay or clearance by the blood.
For a point source the radial dependence of the number densities and alpha
dose components comprises terms proportional to $e^{-r/L_{x}}/(r/L_{x})$. As
noted in (?) typical ranges are $L_{Rn}\sim 0.2-0.4$ mm, $L_{Pb}\sim 0.3-0.7$
mm and $L_{Bi}/L_{Pb}\sim 0.1-0.2$.

Another important parameter introduced in (?) is the 212Pb leakage probability
$P_{leak}(Pb)$, defined as the probability that a 212Pb atom released from the
source is cleared from the tumor by the blood before its decay. The leakage
probability therefore reflects the competition between 212Pb radioactive decay
and clearance through the blood, such that:

$P_{leak}(Pb)=\frac{\alpha_{Pb}}{\lambda_{Pb}+\alpha_{Pb}}$ (12)

The typical range of values is $P_{leak}(Pb)\sim 0.2-0.8$.

##  3 Asymptotic and approximate time-dependent solutions for an infinitely
long cylindrical source

As discussed in (?), at long times after source insertion into the tumor the
number densities reach an asymptotic form:
$n_{x}^{asy}(r,z,t)=\widetilde{n}_{x}(r,z)e^{-\lambda_{Ra}t}$. For 220Rn this
condition is satisfied within several minutes throughout the tumor, while for
212Pb and 212Bi the asymptotic form is attained within a few days, depending
on the distance from the source.

In A we provide a derivation of the closed-form asymptotic solution of the DL
model equations for an infinitely long cylindrical source. For 220Rn we get:

$n_{Rn}^{asy}(r,t)=A_{Rn}\;K_{0}\left(\frac{r}{L_{Rn}}\right)e^{-\lambda_{Ra}t}$
(13)

where:

$A_{Rn}=\frac{P_{des}(Rn)\left(\Gamma_{Ra}^{src}(0)/l\right)}{2\pi
D_{Rn}\,(R_{0}/L_{Rn})K_{1}\left(R_{0}/L_{Rn}\right)}$ (14)

In these expressions $K_{0}(\xi)$ and $K_{1}(\xi)$ are modified Bessel
functions of the second kind:

$\displaystyle K_{0}\left(\xi\right)=\int_{0}^{\infty}\frac{\cos\left(\xi
t\right)}{\sqrt{t^{2}+1}}dt$ (15) $\displaystyle
K_{1}(\xi)=-\frac{dK_{0}}{d\xi}$ (16)

$K_{0}(r/L)$ is a steeply-falling function, and is the cylindrical analogue to
$\exp(-r/L)/(r/L)$ appearing in expressions for the number densities and dose
of the point source (?). As shown in A, for 212Pb we have:

$\displaystyle
n_{Pb}^{asy}(r,t)=\left(A_{Pb}\;K_{0}\left(\frac{r}{L_{Rn}}\right)+B_{Pb}\;K_{0}\left(\frac{r}{L_{Pb}}\right)\right)e^{-\lambda_{Ra}t}$
(17)

where:

$\displaystyle\qquad\qquad
A_{Pb}=\frac{L_{Rn}^{2}L_{Pb}^{2}}{L_{Rn}^{2}-L_{Pb}^{2}}\frac{\lambda_{Rn}}{D_{Pb}}A_{Rn}$
(18) $\displaystyle\qquad\qquad
B_{Pb}=\frac{\left(P_{des}^{eff}(Pb)-P_{des}(Rn)\right)\left(\Gamma_{Ra}^{src}(0)/l\right)}{2\pi
D_{Pb}\,(R_{0}/L_{Pb})K_{1}(R_{0}/L_{Pb})}-A_{Pb}\frac{(R_{0}/L_{Rn})K_{1}(R_{0}/L_{Rn})}{(R_{0}/L_{Pb})K_{1}(R_{0}/L_{Pb})}$
(19)

Finally, for 212Bi we get:

$\qquad\qquad
n_{Bi}^{asy}(r,t)=\left(A_{Bi}\;K_{0}\left(\frac{r}{L_{Rn}}\right)+B_{Bi}\;K_{0}\left(\frac{r}{L_{Pb}}\right)+C_{Bi}\;K_{0}\left(\frac{r}{L_{Bi}}\right)\right)e^{-\lambda_{Ra}t}$
(20)

where:

$A_{Bi}=\frac{L_{Rn}^{2}L_{Bi}^{2}}{L_{Rn}^{2}-L_{Bi}^{2}}\frac{\lambda_{Pb}}{D_{Bi}}A_{Pb}$
(21)

$B_{Bi}=\frac{L_{Pb}^{2}L_{Bi}^{2}}{L_{Pb}^{2}-L_{Bi}^{2}}\frac{\lambda_{Pb}}{D_{Bi}}B_{Pb}$
(22)

$C_{Bi}=-\frac{(R_{0}/L_{Rn})K_{1}(R_{0}/L_{Rn})A_{Bi}+(R_{0}/L_{Pb})K_{1}(R_{0}/L_{Pb})B_{Bi}}{(R_{0}/L_{Bi})K_{1}(R_{0}/L_{Bi})}$
(23)

As shown in A the expressions above can also describe the limit of an infinite
line source in the limit $R_{0}/L_{x}\rightarrow 0$.

To approximately account for the buildup stage of the solution, one can assume
that it is uniform throughout the tumor, i.e., independent of the distance
from the source. Under this “0D” temporal approximation, discussed in (?) for
a point source and adapted here for the cylindrical case, one can write:

$n_{Rn}^{0D}(r,t)=A_{Rn}\;K_{0}\left(\frac{r}{L_{Rn}}\right)\left(e^{-\lambda_{Ra}t}-e^{-\lambda_{Rn}t}\right)$
(24)

and

$\displaystyle\qquad n_{Bi}^{0D}(r,t)=$
$\displaystyle\qquad\left(A_{Bi}\;K_{0}\left(\frac{r}{L_{Rn}}\right)+B_{Bi}\;K_{0}\left(\frac{r}{L_{Pb}}\right)+C_{Bi}\;K_{0}\left(\frac{r}{L_{Bi}}\right)\right)\left(e^{-\lambda_{Ra}t}-e^{-(\lambda_{Pb}+\alpha_{Pb})t}\right)$
(25)

Under this approximation, the asymptotic alpha dose components are:

$\qquad
Dose_{\alpha}^{asy}(RnPo;\,r)=\frac{E_{\alpha}(RnPo)}{\rho}\lambda_{Rn}A_{Rn}\;K_{0}\left(\frac{r}{L_{Rn}}\right)\left(\tau_{Ra}-\tau_{Rn}\right)$
(26)

$\displaystyle\qquad Dose_{\alpha}^{asy}(BiPo;\,r)=$
$\displaystyle\qquad\frac{E_{\alpha}(BiPo)}{\rho}\lambda_{Bi}\left(A_{Bi}\;K_{0}\left(\frac{r}{L_{Rn}}\right)+B_{Bi}\;K_{0}\left(\frac{r}{L_{Pb}}\right)+C_{Bi}\;K_{0}\left(\frac{r}{L_{Bi}}\right)\right)\left(\tau_{Ra}-\tau_{Pb}^{eff}\right)$

where $\tau_{Ra}=1/\lambda_{Ra}$, $\tau_{Rn}=1/\lambda_{Rn}$ and
$\tau_{Pb}^{eff}=1/(\lambda_{Pb}+\alpha_{Pb})$. The error introduced by the 0D
approximation is of the order of the ratio between mean lifetimes of 220Rn and
212Pb and that of 224Ra, i.e. $\tau_{Rn}/\tau_{Ra}\sim 10^{-4}$ and
$\tau_{Pb}/\tau_{Ra}\sim 0.1$, respectively.

##  4 Finite-element one-dimensional time-dependent solution for an infinite
cylindrical source

###  4.1 Numerical scheme

A complete time-dependent solution to the DL model can be done numerically
using a finite-element approach. For the one-dimensional case, i.e., infinite
cylindrical or line source along the $z$ axis, we solve equations (1-3)
setting $\partial^{2}n_{x}/\partial z^{2}=0$. The solution therefore depends
solely on the radial coordinate $r$. We divide our domain into concentric
cylindrical shells, enumerated $i=1...N_{r}$, of equal radial width $\Delta
r$. The radius of the source is $R_{0}$; we choose $\Delta r$ such that
$R_{0}/\Delta r$ is an integer number, and $\Delta r$ is considerably smaller
than $L_{Rn}$ and $L_{Pb}$ ($L_{Bi}$ has a negligible effect on the solution
and therefore should not constrain $\Delta r$). The central radius of the
$i$-th shell is:

$r_{i}=R_{0}+(i-\frac{1}{2})\Delta r$ (28)

We employ a fully implicit scheme (?), thereby assuring the solution is
stable. The time steps are changed adaptively according to the relative change
in the solution between the current step and the previous one, as explained
below.

The average number densities in the $i$-th shell are $n_{Rn,i}$, $n_{Pb,i}$
and $n_{Bi,i}$. We enumerate the time steps by $p$. For shells away from the
source surface, with $1<i\leq N_{r}$, the DL model equations take the discrete
implicit form:

$\displaystyle\qquad\frac{n^{(p+1)}_{x,i}-n^{(p)}_{x,i}}{\Delta t}=$
$\displaystyle\qquad
D_{x}\bigg{(}\frac{n^{(p+1)}_{x,i+1}+n^{(p+1)}_{x,i-1}-2n^{(p+1)}_{x,i}}{\Delta
r^{2}}+\frac{1}{r_{i}}\frac{n^{(p+1)}_{x,i+1}-n^{(p+1)}_{x,i-1}}{2\Delta
r}\bigg{)}-(\lambda_{x}+\alpha_{x})n^{(p+1)}_{x,i}+s^{(p+1)}_{x,i}$

where $x$ stands for Rn, Pb and Bi and $\alpha_{Rn}=0$. Outside our domain we
set the number densities to zero, such that in eq. (4.1) for $i=N_{r}$
$n_{x,i+1}=0$.

As shown in B, for the $i=1$ shell, immediately outside of the source, using
the boundary conditions, eq. (4)-(6), gives:

$\displaystyle\qquad\frac{n^{(p+1)}_{x,1}-n^{(p)}_{x,1}}{\Delta
t}=\frac{D_{x}}{\Delta r^{2}}\left(\frac{1+\Delta r/R_{0}}{1+\Delta
r/2R_{0}}\right)\left(n^{(p+1)}_{x,2}-n^{(p+1)}_{x,1}\right)-(\lambda_{x}+\alpha_{x})n^{(p+1)}_{x,1}+s^{(p+1)}_{x,1}$

The source terms $s^{p+1}_{x,i}$ appearing in eq. (4.1) and (4.1) are:

$\displaystyle\qquad
s_{Rn,i}^{p+1}=\frac{P_{des}(Rn)\left(\Gamma_{Ra}^{src}(0)/l\right)e^{-\lambda_{Ra}t_{p+1}}}{2\pi
R_{0}\Delta r\left(1+\Delta r/2R_{0}\right)}\delta_{i,1}$ (31)
$\displaystyle\qquad
s_{Pb,i}^{p+1}=\frac{\left(P_{des}^{eff}(Pb)-P_{des}(Rn)\right)\left(\Gamma_{Ra}^{src}(0)/l\right)e^{-\lambda_{Ra}t_{p+1}}}{2\pi
R_{0}\Delta r\left(1+\Delta
r/2R_{0}\right)}\delta_{i,1}+\lambda_{Rn}n_{Rn,i}^{p+1}$ (32)
$\displaystyle\qquad s_{Bi,i}^{p+1}=\lambda_{Pb}n_{Pb,i}^{p+1}$ (33)

where $\delta_{i,1}=1$ for $i=1$ and zero otherwise. Rearranging eq. (4.1) and
(4.1), we get the general form:

$n^{(p)}_{x,i}+s^{(p+1)}_{x,i}\Delta
t=A_{i,i-1}^{(x)}n^{(p+1)}_{x,i-1}+A_{i,i}^{(x)}n^{(p+1)}_{x,i}+A_{i,i+1}^{(x)}n^{(p+1)}_{x,i+1}$
(34)

The matrix coefficients introduced in eq. (34) depend on the value of $i$,
reflecting the boundary conditions for $i=1$ and $i=N_{r}$. Retaining terms up
to first order in $\Delta r/r_{i}$ the different cases are summarized below:

$\displaystyle\qquad\qquad A_{i,i-1}^{(x)}$ $\displaystyle=-\frac{D_{x}\Delta
t}{\Delta r^{2}}\left(1-\frac{\Delta r}{2r_{i}}\right)$ $\displaystyle 1<i\leq
N_{r}$ (35) $\displaystyle\qquad\qquad A_{i,i}^{(x)}$
$\displaystyle=1+\frac{D_{x}\Delta t}{\Delta r^{2}}\left(1+\frac{\Delta
r}{2r_{i}}\right)+\left(\lambda_{x}+\alpha_{x}\right)\Delta t\qquad$
$\displaystyle i=1$ (36) $\displaystyle=1+\frac{2D_{x}\Delta t}{\Delta
r^{2}}+\left(\lambda_{x}+\alpha_{x}\right)\Delta t$ $\displaystyle 1<i\leq
N_{r}$ $\displaystyle\qquad\qquad A_{i,i+1}^{(x)}$
$\displaystyle=-\frac{D_{x}\Delta t}{\Delta r^{2}}\left(1+\frac{\Delta
r}{2r_{i}}\right)$ $\displaystyle 1\leq i<N_{r}$ (37)

with $r_{i}$ given in eq. (28). Writing eq. (34) in matrix form:

which can be written as $\mathbf{n}^{(p)}_{x}+\mathbf{s}^{(p+1)}_{x}\Delta
t=\mathbf{A}_{x}\mathbf{n}^{(p+1)}_{x}$. Multiplying on the left by the
inverse of $\mathbf{A}_{x}$, we get:

$\mathbf{n}^{(p+1)}_{x}=\mathbf{A}_{x}^{-1}(\mathbf{n}^{(p)}_{x}+\mathbf{s}^{(p+1)}_{x}\Delta
t),$ (39)

which completes the solution for the $p+1$ step. Note that although the source
terms are calculated for the $p+1$ step they are, in fact, known when the
matrix equations are solved. The reason is that in the $p+1$ step we first
solve for 220Rn, then for 212Pb and lastly for 212Bi. The source term for
220Rn depends only on time, those of 212Pb are found using the $p+1$ solution
for 220Rn, and those of 212Bi – using the $p+1$ solution for 212Pb. Another
point to take into account is that since $\Delta t$ is changed along the
calculation and the matrix coefficients depend on $\Delta t$, they must be
updated accordingly in each step.

The alpha dose components are also updated in each step:

$\qquad\qquad
Dose_{\alpha}^{(p+1)}(RnPo;\,i)=Dose_{\alpha}^{(p)}(RnPo;\,i)+\frac{E_{\alpha}(RnPo)}{\rho}\lambda_{Rn}n_{Rn,i}^{(p+1)}\Delta
t$ (40)

$\qquad\qquad
Dose_{\alpha}^{(p+1)}(BiPo;\,i)=Dose_{\alpha}^{(p)}(BiPo;\,i)+\frac{E_{\alpha}(BiPo)}{\rho}\lambda_{Bi}n_{Bi,i}^{(p+1)}\Delta
t$ (41)

At the end of the $p+1$ step, $\Delta t$ is updated based on the relative
change in the solution. This can be done in a number of ways. A particular
choice, implemented here, was to consider the relative change in the total
dose (sum of the $RnPo$ and $BiPo$ contributions) at a particular point of
interest $r_{i_{0}}$ (e.g., at 2 mm):

$\qquad\qquad\Delta t_{new}={\Delta
t}\cdot\frac{\epsilon_{tol}}{\left(Dose_{\alpha}^{(p+1)}(tot;i_{0})-Dose_{\alpha}^{(p)}(tot;i_{0})\right)/Dose_{\alpha}^{(p)}(tot;i_{0})}$
(42)

where $\epsilon_{tol}$ is a preset tolerance parameter. For practicality, one
can further set upper and lower limits on $\Delta t$ to balance between
calculation time and accuracy. Although the initial time step should be small
compared to 220Rn half-life, its particular value has little effect on the
accuracy of the calculated asymptotic dose.

###  4.2 DART1D: tests and comparisons

The finite-element scheme described above was implemented in MATLAB in a code
named “DART1D”. The solution of eq. (39), the critical part of the
calculation, was done using a tridiagonal matrix solver employing the Thomas
algorithm (?). The use of this solver was found to be $\sim 4$ times faster
than MATLAB’s mldivide (‘\’) tool, which was, in turn, about 3-fold faster
than inverting the matrix using inv(A). The code was found to converge to sub-
percent level for a modest choice of the discretization parameter values. For
example, setting $\epsilon_{tol}=10^{-2}$, $\Delta r=0.02$ mm, and $\Delta
t_{0}=0.1$ s (for a domain radius $R_{max}=7$ mm and a treatment duration of
30 d) resulted, with a run-time of $\sim 0.5$ s, in doses which were $\sim
0.5\%$ away from those obtained with $\epsilon_{tol}=10^{-4}$, $\Delta r=0.01$
mm and $\Delta t_{0}=0.1$ s, with a run-time of $\sim 3$ min (both on a modern
laptop computer with an Intel i7 processor and 16 GB RAM memory). The latter,
more accurate run, was within $7\cdot 10^{-4}$ of the 0D approximation for the
220Rn+216Po dose.

Figures 2 and 3 examine several aspects of the numerical solution. Figure 2
shows the DART1D solution in comparison with the asymptotic expressions eq.
(13) and (17). On the left we show the DART1D time-dependent 212Pb number
density at a distance of 2 mm from the source, along with the corresponding
asymptotic solution. On the right, we show the ratio between the numerical and
asymptotic solutions for 220Rn, $f_{Rn}\equiv n_{Rn}/n_{Rn}^{asy}$, plotted
for varying distances from the source. The distances are given in units of the
220Rn diffusion length, $r^{*}\equiv r/L_{Rn}$, and the time in units of
$1/(\lambda_{Rn}-\lambda_{Ra})$ which is roughly the mean 220Rn lifetime,
$t^{*}\equiv(\lambda_{Rn}-\lambda_{Ra})t$. The numerical solutions converge to
the asymptotic ones with a delay that increases with the distance from the
source. For 220Rn this occurs on the scale of minutes, while for 212Pb – over
a few days. The adaptive time step allows DART1D to handle both transients
efficiently.

Figure 3A shows the DART1D 220Rn+216Po and 212Bi/212Po alpha dose components
calculated for the case $L_{Rn}=0.3$ mm, $L_{Pb}=0.6$ mm, $L_{Bi}=0.1L_{Pb}$,
$\alpha_{Pb}=\lambda_{Pb}$ (i.e., $P_{leak}(Pb)=0.5$), and $\alpha_{Bi}=0$.
The source radius is $R_{0}=0.35$ mm, the 224Ra activity is 3 $\mu$Ci/cm and
the desorption probabilities are $P_{des}(Rn)=0.45$ and
$P_{des}^{eff}(Pb)=0.55$. The dose components are given at $t=30$ d post
treatment. The numerical calculation is compared to the 0D approximations, eq.
(26) and (3). The assumption of zero number density outside the calculation
domain results in a departure from the expected solution about two diffusion
lengths away from the boundary: $\sim 0.5$ mm for 220Rn and $\sim 1$ mm for
212Bi and 212Po, whose spatial distribution is governed by the 212Pb diffusion
length. This indicates that the radial extent of the calculation domain should
be about 10 times larger than the largest diffusion length of the problem.
Figure 3B shows the ratio between the DART1D-calculated dose components and
the corresponding 0D approximations. Except for the edge effect at
$r\rightarrow R_{max}$, the numerical solution for 220Rn+216Po coincides with
the 0D approximation to better than $1\cdot 10^{-3}$ for
$\epsilon_{tol}=10^{-4}$, $\Delta r=0.01$ mm and $\Delta t_{0}=0.1$ s up to
$r\sim 5$ mm. For 212Bi/212Po the 0D approximation underestimates the dose at
$r<1$ mm and overestimates it at larger distances because of the increasing
delay in buildup of 212Pb as a function of $r$. The error is $\sim 5-10\%$ at
therapeutically relevant distances from the source (around 2-3 mm).

Figure 2: Comparison between the DART1D and asymptotic solutions for the 212Pb
number density 2 mm from the source (A) and the ratio between the DART1D and
asymptotic solutions of the 220Rn number density at various distances from the
source axis (B). The distance and time are normalized to 220Rn diffusion
length and mean lifetime: $r^{*}\equiv r/L_{Rn}$,
$t^{*}\equiv(\lambda_{Rn}-\lambda_{Ra})t$.

Figure 3: Comparison between the DART1D 220Rn+216Po and 212Bi/212Po alpha
doses and the 0D approximation for an infinite cylindrical source. (A) Dose
values, (B) DART1D/0D approximation ratios.

##  5 Finite-element two-dimensional time-dependent solution for a finite
cylindrical source

###  5.1 Numerical scheme

We move now to two dimensions to treat a cylindrical source (“seed”) of radius
$R_{0}$ and finite length $l$. The source lies along the $z$ axis with $z=0$
at its mid plane. We solve the DL equations over a cylindrical domain
extending from $r=0$ to $r=R_{max}$ and from $z=-Z_{max}$ to $z=+Z_{max}$.
Both $R_{max}$ and $Z_{max}-l/2$ should be much larger than the largest
diffusion length of the problem. The domain comprises ring elements of equal
radial width $\Delta r$ and equal $z$-width $\Delta z$. We choose $\Delta r$
and $\Delta z$ such that $R_{0}/\Delta r$ and $l/(2\Delta z)$ are integer
numbers, with $\Delta r$ and $\Delta z$ much smaller than $L_{Rn}$ and
$L_{Pb}$. We enumerate the rings by $i,j$, where $i=1...N_{r}$ and
$j=1...N_{z}$. Elements with $i=1$ are on-axis, while $j=1$ at the bottom of
the cylindrical domain, and $j=N_{z}$ at the top. Unlike the 1D case, where
the source is infinitely long and we only consider points with $r>R_{0}$, for
a finite seed in 2D we must also solve the equations for points above and
below the seed, with $r<R_{0}$ and $|z|>\frac{1}{2}l$. As for the 1D case, the
radius and $z$-coordinate of the $i,j$ ring, $r_{i},z_{j}$, are defined at the
center of its $rz$ cross section. For the innermost $i=1$ rings
$r_{1}=\frac{1}{2}\Delta r$. Points inside the seed, i.e., in rings with
$r_{i}\leq R_{0}-\Delta r/2$, and $|z_{j}|\leq\frac{1}{2}l-\frac{1}{2}\Delta
z$, have zero number densities of 220Rn, 212Pb and 212Bi.

Discretization of eq. (1-3) in 2D yields, for interior ring elements in the
cylindrical domain (outside of the seed and not touching its wall or bases,
and with $i>1$):

$\eqalign{\frac{n^{(p+1)}_{x,i,j}-n^{(p)}_{x,i,j}}{\Delta t}=\cr
D_{x}\bigg{(}\frac{n^{(p+1)}_{x,i+1,j}+n^{(p+1)}_{x,i-1,j}-2n^{(p+1)}_{x,i,j}}{\Delta
r^{2}}+\frac{1}{r_{i}}\frac{n^{(p+1)}_{x,i+1,j}-n^{(p+1)}_{x,i-1,j}}{2\Delta
r}\cr+\frac{n^{(p+1)}_{x,i,j+1}+n^{(p+1)}_{x,i,j-1}-2n^{(p+1)}_{x,i,j}}{\Delta
z^{2}}\bigg{)}-(\lambda_{x}+\alpha_{x})n^{(p+1)}_{x,i,j}+s_{x,i,j}^{(p+1)}}$
(43)

Eq. (43) holds also for ring elements on the external surfaces of the domain,
with $i=N_{r}$, $j=1$ or $j=N_{z}$, as the number densities for rings with
$i=N_{r}+1$, $j=0$ or $j=N_{z}+1$ are all zero. For ring elements on-axis
($i=1$), above or below the seed, we require $\left(\partial n_{x}/\partial
r\right)_{r=0}=0$. Since the 224Ra activity is confined to the seed wall, the
$z$-component of the current density is set to zero on the seed bases, i.e.,
$\left(\partial n_{x}\partial z\right)_{z=\pm l/2}=0$. Defining $i_{s}$ as the
radial index of ring elements touching the seed wall (i.e.,
$r_{i_{s}}=R_{0}+\Delta r/2$), for ring elements with $|z_{j}|\leq l/2-\Delta
z/2$, eq. (43) becomes, similarly to the 1D case:

$\displaystyle\qquad\frac{n^{(p+1)}_{x,i_{s},j}-n^{(p)}_{x,i_{s},j}}{\Delta
t}=$ $\displaystyle\frac{D_{x}}{\Delta r^{2}}\left(\frac{1+\Delta
r/R_{0}}{1+\Delta
r/2R_{0}}\right)\left(n^{(p+1)}_{x,i_{s}+1,j}-n^{(p+1)}_{x,i_{s},j}\right)+$
$\displaystyle\frac{D_{x}}{\Delta
z^{2}}\left(n^{(p+1)}_{x,i_{s},j+1}+n^{(p+1)}_{x,i_{s},j-1}-2n^{(p+1)}_{x,i_{s},j}\right)-(\lambda_{x}+\alpha_{x})n^{(p+1)}_{x,i_{s},j}+s^{(p+1)}_{x,i_{s},j}$

The source terms in eq. (43) and (5.1) are similar to the 1D case, with the
additional requirement that $|z_{j}|<l/2$:

$\qquad
s_{Rn,i,j}^{p+1}=\frac{P_{des}(Rn)\left(\Gamma_{Ra}^{src}(0)/l\right)e^{-\lambda_{Ra}t_{p+1}}}{2\pi
R_{0}\Delta r\left(1+\Delta
r/2R_{0}\right)}\delta_{i,i_{s}}\cdot\left(\frac{1-\mathrm{sign}(|z_{j}|-l/2)}{2}\right)$
(45)

$\displaystyle\qquad
s_{Pb,i,j}^{p+1}=\frac{\left(P_{des}^{eff}(Pb)-P_{des}(Rn)\right)\left(\Gamma_{Ra}^{src}(0)/l\right)e^{-\lambda_{Ra}t_{p+1}}}{2\pi
R_{0}\Delta r\left(1+\Delta
r/2R_{0}\right)}\delta_{i,i_{s}}\cdot\left(\frac{1-\mathrm{sign}(|z_{j}|-l/2)}{2}\right)$
$\displaystyle+\lambda_{Rn}n_{Rn,i,j}^{p+1}$ (46)

$\qquad s_{Bi,i,j}^{p+1}=\lambda_{Pb}n_{Pb,i,j}^{p+1}$ (47)

In order to solve eq. (43) in matrix form we use linear indexing. We rearrange
the 2D elements $n_{x,i,j}^{(p)}$ and $s_{x,i,j}^{(p)}$ in two column vectors
$\mathbf{\widetilde{n}}_{x}^{(p)}$ and $\mathbf{\widetilde{s}}_{x}^{(p)}$ in
sequential order. We define:

$k(i,j)=(j-1)N_{r}+i$ (48)

$\widetilde{n}_{x,k}^{(p)}=n_{x,i,j}^{(p)}$ (49)

$\widetilde{s}_{x,k}^{(p)}=s_{x,i,j}^{(p)}$ (50)

with $k=1...N_{r}Nz$. Noting that $n_{x,i\pm 1,j}^{(p)}=\widetilde{n}_{x,k\pm
1}^{(p)}$ and $n_{x,i,j\pm 1}^{(p)}=\widetilde{n}_{x,k\pm N_{r}}^{(p)}$, eq.
(43) can be rearranged as:

$\displaystyle\quad\widetilde{n}^{(p)}_{x,k}+\widetilde{s}_{x,k}^{(p+1)}\Delta
t=$ $\displaystyle\quad
M_{k,k-N_{r}}^{(x)}\widetilde{n}^{(p+1)}_{x,k-N_{r}}+M_{k,k-1}^{(x)}\widetilde{n}^{(p+1)}_{x,k-1}+M_{k,k}^{(x)}\widetilde{n}^{(p+1)}_{x,k}+M_{k,k+1}^{(x)}\widetilde{n}^{(p+1)}_{x,k+1}+M_{k,k+N_{r}}^{(x)}\widetilde{n}^{(p+1)}_{x,k+N_{r}}$
(51)

As for the 1D case, the matrix elements appearing in eq. (5.1) depend on the
values of $i,j$ (and therefore $k$), in a manner that satisfies the boundary
conditions. For compactness, we define the following intermediate quantities:

$\displaystyle K_{z}^{(x)}$ $\displaystyle=\frac{D_{x}\Delta t}{\Delta z^{2}}$
(52) $\displaystyle K_{r}^{(x)}$ $\displaystyle=\frac{D_{x}\Delta t}{\Delta
r^{2}}$ (53) $\displaystyle K_{r}^{(x+)}$ $\displaystyle=\frac{D_{x}\Delta
t}{\Delta r^{2}}\left(1+\frac{\Delta r}{2r_{i}}\right)$ (54) $\displaystyle
K_{r}^{(x-)}$ $\displaystyle=\frac{D_{x}\Delta t}{\Delta
r^{2}}\left(1-\frac{\Delta r}{2r_{i}}\right)$ (55) $\displaystyle S_{1}^{(x)}$
$\displaystyle=1+2K_{r}^{(x)}+K_{z}^{(x)}+(\lambda_{x}+\alpha_{x})\Delta t$
(56) $\displaystyle S_{2}^{(x)}$
$\displaystyle=1+2K_{r}^{(x)}+2K_{z}^{(x)}+(\lambda_{x}+\alpha_{x})\Delta t$
(57) $\displaystyle S_{+}^{(x)}$
$\displaystyle=1+K_{r}^{(x+)}+2K_{z}^{(x)}+(\lambda_{x}+\alpha_{x})\Delta t$
(58)

Table 1 lists the expressions for the matrix elements $M_{k,l}^{(x)}$ for all
possible cases for $r_{i}$ and $z_{j}$. With these, we can write eq. (5.1) in
matrix form (with $K\equiv N_{r}N_{z}$):

Table 1: Matrix elements in 2D. Case | $M_{k,k-N_{r}}^{(x)}$ | $M_{k,k-1}^{(x)}$ | $M_{k,k}^{(x)}$ | $M_{k,k+1}^{(x)}$ | $M_{k,k+N_{r}}^{(x)}$
---|---|---|---|---|---
$R_{0}+\frac{\Delta r}{2}<r_{i}<R_{max}-\frac{\Delta r}{2}$ |  |  |  |  |
& $|z_{j}|<Z_{max}-\frac{\Delta z}{2}$ | $-K_{z}^{(x)}$ | $-K_{r}^{(x-)}$ | $S_{2}^{(x)}$ | $-K_{r}^{(x+)}$ | $-K_{z}^{(x)}$
$\frac{\Delta r}{2}<r_{i}\leq R_{0}+\frac{\Delta r}{2}$ |  |  |  |  |
& $\frac{l}{2}+\frac{\Delta z}{2}<|z_{j}|<Z_{max}-\frac{\Delta z}{2}$ | $-K_{z}^{(x)}$ | $-K_{r}^{(x-)}$ | $S_{2}^{(x)}$ | $-K_{r}^{(x+)}$ | $-K_{z}^{(x)}$
$r_{i}=R_{0}+\frac{\Delta r}{2}$ & $|z_{j}|=\frac{l}{2}+\frac{\Delta z}{2}$ | $-K_{z}^{(x)}$ | $-K_{r}^{(x-)}$ | $S_{2}^{(x)}$ | $-K_{r}^{(x+)}$ | $-K_{z}^{(x)}$
$r_{i}=R_{0}+\frac{\Delta r}{2}$ & $|z_{j}|<\frac{l}{2}$ | $-K_{z}^{(x)}$ | 0 | $S_{+}^{(x)}$ | $-K_{r}^{(x+)}$ | $-K_{z}^{(x)}$
$\frac{\Delta r}{2}<r_{i}<R_{0}$ & $z_{j}=\frac{l}{2}+\frac{\Delta z}{2}$ | 0 | $-K_{r}^{(x-)}$ | $S_{1}^{(x)}$ | $-K_{r}^{(x+)}$ | $-K_{z}^{(x)}$
$\frac{\Delta r}{2}<r_{i}<R_{0}$ & $z_{j}=-\frac{l}{2}-\frac{\Delta z}{2}$ | $-K_{z}^{(x)}$ | $-K_{r}^{(x-)}$ | $S_{1}^{(x)}$ | $-K_{r}^{(x+)}$ | 0
$r_{i}=\frac{\Delta r}{2}$ |  |  |  |  |
& $\frac{l}{2}+\frac{\Delta z}{2}<|z_{j}|<Z_{max}-\frac{\Delta z}{2}$ | $-K_{z}^{(x)}$ | 0 | $S_{2}^{(x)}$ | $-K_{r}^{(x+)}$ | $-K_{z}^{(x)}$
$r_{i}=R_{max}-\frac{\Delta r}{2}$ |  |  |  |  |
& $|z_{j}|<Z_{max}-\frac{\Delta z}{2}$ | $-K_{z}^{(x)}$ | $-K_{r}^{(x-)}$ | $S_{2}^{(x)}$ | 0 | $-K_{z}^{(x)}$
$\frac{\Delta r}{2}<r_{i}<R_{max}-\frac{\Delta r}{2}$ |  |  |  |  |
& $z_{j}=Z_{max}-\frac{\Delta z}{2}$ | $-K_{z}^{(x)}$ | $-K_{r}^{(x-)}$ | $S_{2}^{(x)}$ | $-K_{r}^{(x+)}$ | 0
$\frac{\Delta r}{2}<r_{i}<R_{max}-\frac{\Delta r}{2}$ |  |  |  |  |
& $z_{j}=-Z_{max}+\frac{\Delta z}{2}$ | 0 | $-K_{r}^{(x-)}$ | $S_{2}^{(x)}$ | $-K_{r}^{(x+)}$ | $-K_{z}^{(x)}$
$r_{i}=\frac{\Delta r}{2}\;\&\;z_{j}=Z_{max}-\frac{\Delta z}{2}$ | $-K_{z}^{(x)}$ | 0 | $S_{2}^{(x)}$ | $-K_{r}^{(x+)}$ | 0
$r_{i}=\frac{\Delta r}{2}\;\&\;z_{j}=-Z_{max}+\frac{\Delta z}{2}$ | 0 | 0 | $S_{2}^{(x)}$ | $-K_{r}^{(x+)}$ | $-K_{z}^{(x)}$
$r_{i}=\frac{\Delta r}{2}\;\&\;z_{j}=\frac{l}{2}+\frac{\Delta z}{2}$ | 0 | 0 | $S_{2}^{(x)}$ | $-K_{r}^{(x+)}$ | $-K_{z}^{(x)}$
$r_{i}=\frac{\Delta r}{2}\;\&\;z_{j}=-\frac{l}{2}-\frac{\Delta z}{2}$ | $-K_{z}^{(x)}$ | 0 | $S_{2}^{(x)}$ | $-K_{r}^{(x+)}$ | 0
$r_{i}=R_{max}-\frac{\Delta r}{2}$ & $z_{j}=Z_{max}-\frac{\Delta z}{2}$ | $-K_{z}^{(x)}$ | $-K_{r}^{(x-)}$ | $S_{2}^{(x)}$ | 0 | 0
$r_{i}=R_{max}-\frac{\Delta r}{2}$ & $z_{j}=-Z_{max}+\frac{\Delta z}{2}$ | 0 | $-K_{r}^{(x-)}$ | $S_{2}^{(x)}$ | 0 | $-K_{z}^{(x)}$
$r_{i}<R_{0}$ & $|z_{j}|<\frac{l}{2}$ | 0 | 0 | 1 | 0 | 0

or, equivalently:
$\mathbf{\widetilde{n}}^{(p)}_{x}+\mathbf{\widetilde{s}}^{(p+1)}_{x}\Delta
t=\mathbf{M}_{x}\mathbf{\widetilde{n}}^{(p+1)}_{x}$. As for the 1D case, we
multiply on the left by the inverse of $\mathbf{M}_{x}$, getting
$\mathbf{\widetilde{n}}^{(p+1)}_{x}=\mathbf{M}_{x}^{-1}(\mathbf{\widetilde{n}}^{(p)}_{x}+\mathbf{\widetilde{s}}^{(p+1)}_{x}\Delta
t)$. We now run over all possible values of $i,j$ and update
$n_{x,i,j}^{(p+1)}=\widetilde{n}_{x,k}^{(p+1)}$. Once the new number densities
are known in all ring elements, we update the alpha dose components:

$\qquad\qquad
Dose_{\alpha}^{(p+1)}(RnPo;\,i,j)=Dose_{\alpha}^{(p)}(RnPo;\,i,j)+\frac{E_{\alpha}(RnPo)}{\rho}\lambda_{Rn}n_{Rn,i,j}^{(p+1)}\Delta
t$ (60)

$\qquad\qquad
Dose_{\alpha}^{(p+1)}(BiPo;\,i,j)=Dose_{\alpha}^{(p)}(BiPo;\,i,j)+\frac{E_{\alpha}(BiPo)}{\rho}\lambda_{Bi}n_{Bi,i,j}^{(p+1)}\Delta
t$ (61)

As for the 1D case, the time step can be modified in many ways. Here we chose
to update it according to the relative change in the overall activity (sum
over all isotopes in all ring elements).

###  5.2 DART2D: tests and comparisons

The 2D numerical scheme described above was implemented in MATLAB in a code
named “DART2D”. The code takes roughly 0.5 h to run on a modern laptop (Intel
i7 processor with 16 GB RAM) for $\Delta r=0.005$ mm, $\Delta z=0.05$ mm,
$\epsilon_{tol}=0.01$, $\Delta t_{0}=0.1$ s, $R_{max}=7$ mm, $Z_{max}=10$ mm
and a treatment time of 40 d. The most demanding process is the calculation of
$\mathbf{M}^{-1}$. Since $\mathbf{M}$ is a sparse diagonal matrix, we used
MATLAB’s spdiags() function, which reduces memory requirements by saving only
the diagonal non-zero elements of $\mathbf{M}$, and allows the code to run
more efficiently.

Figure 4 shows the dependence of the adaptive time step on time, up to $\sim
11$ days. The initial time step is 0.1 s, capturing the 220Rn buildup with
high accuracy (this is also the minimal allowed value for $\Delta t$). It then
gradually increases, following the buildup of 212Pb, eventually reaching its
maximal allowed value (here $\sim 1$ h) in the asymptotic phase driven by
224Ra decay rate.

Figure 4: DART2D: variation of the time step vs. time.

The total alpha dose (sum of the 220Rn+216Po and 212Bi/212Po contributions)
accumulated over 40 days of DaRT treatment by a seed of finite dimensions is
displayed in the $rz$ plane in figure 5A. The seed dimensions are $R_{0}=0.35$
mm and $l=10$ mm. The initial 224Ra activity of the seed is 3 $\mu$Ci, with
$P_{des}(Rn)=0.45$ and $P_{des}^{eff}(Pb)=0.55$. The other model parameters
are: $L_{Pb}=0.6$ mm, $L_{Rn}=0.3$ mm, $L_{Bi}=0.06$ mm, $P_{leak}(Pb)=0.5$,
$\alpha_{Bi}=0$. Note that the radial dose profile is nearly unchanged up to
$\sim 1.5$ mm from the seed end. Figure 5B shows the dose profiles along $r$
in the seed mid plane and along $z$ parallel to the seed axis, both set such
that ‘0’ is the seed edge. The dose along the seed axis is smaller by $\sim
30\%$ near the seed edge, with the difference increasing to a factor of $\sim
3$ at 3 mm, compared to that in the mid plane - an important point to consider
in treatment planning. Although a similar effect is observed when
approximating the seed to a finite line source comprised of point-like
segments, as done in (?), this approach leads to significant errors in the
dose because it does not consider the finite diameter of the seed, which
“pushes” the radial dose to larger values.

Figure 5: (A) Total alpha dose accumulated over 40 days of treatment by a DaRT
seed with initial activity of $3\;\mu$Ci 224Ra. The seed radius and length are
0.35 mm and 10 mm, respectively. The other model parameters are given in the
text. (B) Total alpha dose as a function of the distance from the seed edge
along $r$ in the mid plane and along $z$ on the seed axis.

We now compare the results of the full 2D calculation with those obtained
using either the 0D analytical approximations, or the full 1D calculation.
Figure 6 shows these comparisons of the dose profile calculated in the seed
mid planes. On the left we display the comparison for a low-diffusion high-
leakage case, with $L_{Rn}=0.3$ mm, $L_{Pb}=0.3$ mm, and $P_{leak}(Pb)=0.8$,
and on the right – for a high-diffusion low-leakage case, with $L_{Rn}=0.3$
mm, $L_{Pb}=0.6$ mm, and $P_{leak}(Pb)=0.3$. In both cases $L_{Bi}=0.1L_{Pb}$
and $\alpha_{Bi}=0$. The curves show the ratios between the full 2D
calculation with DART2D and those obtained by: (1) approximating the seed to a
finite line comprised of point-like segments and using the 0D approximation,
as was done in (?); (2) using the 0D approximation for an infinite cylindrical
source, eq.(26) and (3), and (3) using the full DART1D calculation.
Approximating the seed to a finite line source leads to an underestimation of
the dose by up to $\sim 80\%$ for both the low- and high-diffusion scenarios.
Using the closed-form 0D approximation for a cylindrical source of radius
$R_{0}$ overestimates the dose at 2-3 mm by $\sim 1-2\%$ for the low-
diffusion/high-leakage case and $\sim 5-10\%$ for the high-diffusion/low-
leakage scenario. In contrast, the full numerical solution (DART1D) for a
cylindrical source provides accurate results (on the scale of $0.3\%$) when
compared to the 2D calculation.

Figure 7 shows the dose calculated for a hexagonal seed lattice of parallel
seeds with a grid spacing of 4 mm. As before, the seed radius is 0.35 mm, its
224Ra activity is 3 $\mu$Ci (over 1 cm length), $P_{des}(Rn)=0.45$,
$P_{des}^{eff}(Pb)=0.55$. The calculations are for the low-diffusion high-
leakage and high-diffusion low-leakage cases defined above. The calculation is
done for both the full 2D solution and the 0D line-source approximation. The
dose at the mid point between three adjacent seeds is 75 / 15 Gy for the
accurate 2D calculation (high- / low-diffusion, respectively) and 57 / 8 Gy
for the line source approximation, emphasizing the need to consider the finite
diameter of the seed.

Figure 6: Ratios between the total alpha dose in the seed mid plane calculated
by DART2D and those calculated using the 0D line source approximation, the 0D
infinite cylinder approximation and the full 1D calculation for an infinite
cylindrical source using DART1D. (A) Low-diffusion, high-leakage case; (B)
High-diffusion, low leakage.

Figure 7: Lattice total alpha dose map comparisons. (A) $L_{Rn}=0.3$ mm,
$L_{Pb}=0.3$ mm, $P_{leak}(Pb)=0.8$, line source 0D approximation. (B)
$L_{Rn}=0.3$ mm, $L_{Pb}=0.6$ mm, $P_{leak}(Pb)=0.3$ line source 0D
approximation. (C) Same parameters as A, full 2D calculation with DART2D. (D)
same parameters as B, full 2D calculation with DART2D. In all figures the grid
spacing is 4 mm, seed 224Ra activity 3 $\mu$Ci/cm, $P_{des}(Rn)=0.45$,
$P_{des}^{eff}(Pb)=0.55$. Regardless of the model parameters, taking into
account the finite diameter of the seed with DART2D leads to significantly
higher doses compared to the line source 0D approximation.

##  6 Summary and conclusions

In this work we provided closed-form approximations and numerical finite-
element schemes for calculating the alpha particle dose of DaRT seeds of
finite diameter and length, expanding the discussion of a previous publication
(?) which modeled seeds as idealized line sources comprised of a collection of
point-like segments. The intention of the detailed description of the
numerical schemes was to enable other research groups to develop their own
codes and test them against ours.

The full 2D calculation presented here was shown to better describe the alpha
particle dose field than the point source-based approach, considering the
effect of finite seed diameter and the dose fall-off along the seed axis. We
therefore recommend using it, rather than the simpler approximation, as the
basis for preparing dose lookup tables for DaRT treatment planning.

The intermediate derivation of closed-form approximations and a finite-element
scheme for an infinite cylindrical source can serve two purposes. First, since
the full 1D solution (“DART1D”) coincides to high accuracy with the 2D
calculation on the seed mid plane, it can be used to validate it. In this
respect, it is important to note that due to the rapid buildup of the 220Rn
activity towards its asymptotic form, the accuracy of the closed-form 0D
approximation to the 220Rn+216Po alpha dose is better than $10^{-3}$. It can
therefore be used for validating the calculation of this part of the dose in
both the 1D and 2D finite element schemes.

Second, given that the full 2D calculation is computationally intensive and
can take a long time to complete, the accuracy of DART1D in representing the
dose at the seed mid plane allows using it safely for parameter scans to
understand the properties of DaRT seed lattices. Here it is important to note
that deviations from the mid plane dose profile occur only within $\sim 1.5$
mm from the seed end, as evident from figure 5.

The main aim of the text was to outline the numerical scheme of the 1D and 2D
solutions, rather than develop computer codes with optimized performance. The
most demanding aspect of the calculation is the inversion of the coefficient
matrix, which can be optimized and coded more efficiently than presented here.
For example, when implemented in Fortran, DART1D completes in $\sim$6 s a
calculation which requires $\sim 3$ m to complete in MATLAB.

It is important to emphasize that while the DL model provides a pragmatic
description of the DaRT alpha dose, it is by no means a complete theory. In
particular, it does not take into account the possibility of convective
effects and nonuniformity in real tumors, which comprise both necrotic and
viable regions, with the former evolving as a result of the DaRT treatment
itself. However, in spite of its limitations, the diffusion leakage model can
provide a quantitative guide for selecting a starting point for treatment
planning in clinical trials in terms of seed activity and spacing. This was
demonstrated in the first clinical trial (?), where 2 $\mu$Ci seeds were
inserted at $\sim 5$ mm spacing, based on the DL model prediction that this
would provide a nominal alpha dose of $>10$ Gy throughout the treated volume,
and where $\sim 80\%$ of the treated tumors displayed a complete response.

The predictions of the DL model depend critically on the values of its
parameters, in particular, the diffusion lengths of 220Rn and 212Pb. Extensive
experimental work on measuring them in mouse tumors is ongoing and will be
covered in separate publications.

##  Appendix A Derivation of the asymptotic solutions for infinite cylindrical
and line sources

In this appendix we provide a detailed derivation of the asymptotic solutions
of the diffusion-leakage model equations for the number densities of 220Rn,
212Pb and 212Bi, for an infinite cylindrical source and infinite line source.

###  A.1 220Rn

The 220Rn diffusion equation for the case of an infinite cylindrical source of
radius $R_{0}$ in a homogeneous and isotropic medium, in cylindrical
coordinates, is:

$\frac{\partial n_{Rn}}{\partial t}=\frac{D_{Rn}}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial n_{Rn}}{\partial r}\right)+s_{Rn}-\lambda_{Rn}n_{Rn}$
(62)

Assuming no 224Ra release from the source $s_{Rn}(r,t)=0$, and 220Rn enters
the tumor by direct release from the source surface, with the following
boundary condition at $r\rightarrow R_{0}$:

$\lim_{r\rightarrow R_{0}}2\pi
rj_{Rn}(r,t)=P_{des}(Rn)\frac{\Gamma_{Ra}^{src}(0)}{l}e^{-\lambda_{Ra}t}$ (63)

where $\Gamma_{Ra}^{src}(0)/l$ is the initial 224Ra activity per unit length
of the source, and $j_{Rn}=-D_{Rn}\frac{\partial n_{Rn}}{\partial r}$.
Substituting the asymptotic form
$n_{Rn}^{asy}(r,t)=\widetilde{n}_{Rn}(r)e^{-\lambda_{Ra}t}$ in eq. (62) leads
to:

$\frac{d^{2}\widetilde{n}_{Rn}}{dr^{2}}+\frac{1}{r}\frac{d\widetilde{n}_{Rn}}{dr}-\frac{1}{L_{Rn}^{2}}\widetilde{n}_{Rn}=0$
(64)

where $L_{Rn}$ is defined in eq. (9). Defining $\xi=r/L_{Rn}$ gives:

$\xi^{2}\frac{d^{2}\widetilde{n}_{Rn}}{d\xi^{2}}+\xi\frac{d\widetilde{n}_{Rn}}{d\xi}-\left(\xi^{2}+n^{2}\right)\widetilde{n}_{Rn}=0\qquad
n=0$ (65)

This has the form of a modified Bessel equation (?), for which the solution
is:

$\widetilde{n}_{Rn}(\xi)=A_{Rn}K_{0}(\xi)+B_{Rn}I_{0}(\xi)$ (66)

$I_{0}(\xi)$ is a modified Bessel function of the first kind. Since it
diverges for $\xi\rightarrow\infty$, we must have $B_{Rn}=0$. $K_{0}(\xi)$ is
a modified Bessel function of the second kind, for which:

$K_{0}(\xi)=\int_{0}^{\infty}\frac{\cos\left(\xi t\right)}{\sqrt{t^{2}+1}}dt$
(67)

It vanishes in the limit $\xi\rightarrow\infty$ and hence:

$\widetilde{n}_{Rn}(\xi)=A_{Rn}K_{0}(\xi)$ (68)

The modified Bessel functions of the second kind have the property that:

$\xi\frac{dK_{n}}{d\xi}=nK_{n}(\xi)-\xi K_{n+1}(\xi)$ (69)

Thus:

$\frac{dK_{0}}{d\xi}=-K_{1}(\xi)$ (70)

Using this, the current can be expressed as:

$j_{Rn}(r,t)=\frac{D_{Rn}}{L_{Rn}}A_{Rn}K_{1}\left(\frac{r}{L_{Rn}}\right)e^{-\lambda_{Ra}t}$
(71)

Substituting (71) into the boundary condition (63) yields:

$A_{Rn}=\frac{P_{des}(Rn)\left(\Gamma_{Ra}^{src}(0)/l\right)}{2\pi
D_{Rn}\,(R_{0}/L_{Rn})K_{1}\left(R_{0}/L_{Rn}\right)}$ (72)

Finally:

$n_{Rn}^{asy}(r,t)=\frac{P_{des}(Rn)\left(\Gamma_{Ra}^{src}(0)/l\right)}{2\pi
D_{Rn}\,(R_{0}/L_{Rn})K_{1}\left(R_{0}/L_{Rn}\right)}K_{0}\left(\frac{r}{L_{Rn}}\right)e^{-\lambda_{Ra}t}$
(73)

For $\xi\rightarrow 0$ (?):

$\lim_{\xi\rightarrow 0}\,\xi K_{1}(\xi)=1$ (74)

Therefore, for the case of an ideal line source ($R_{0}/L_{Rn}\rightarrow 0$):

$n_{Rn}^{asy}(r,t)=\frac{P_{des}(Rn)\left(\Gamma_{Ra}^{src}(0)/l\right)}{2\pi
D_{Rn}}K_{0}\left(\frac{r}{L_{Rn}}\right)e^{-\lambda_{Ra}t}$ (75)

###  A.2 212Pb

The diffusion-leakage equation for 212Pb for the infinite cylindrical source
geometry is:

$\frac{\partial n_{Pb}}{\partial t}=\frac{D_{Pb}}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial n_{Pb}}{\partial
r}\right)-\lambda_{Pb}n_{Pb}-\alpha_{Pb}n_{Pb}+\lambda_{Rn}n_{Rn}$ (76)

with the boundary condition:

$\lim_{r\rightarrow R_{0}}2\pi
rj_{Pb}(r,t)=\left(P_{des}^{eff}(Pb)-P_{des}(Rn)\right)\frac{\Gamma_{Ra}^{src}(0)}{l}e^{-\lambda_{Ra}t}$
(77)

where $j_{Pb}=-D_{Pb}\frac{\partial n_{Pb}}{\partial r}$. As before, we are
looking for a solution of the form:

$n_{Pb}^{asy}\left(r,t\right)=\widetilde{n}_{Pb}\left(r\right)e^{-\lambda_{Ra}t}$
(78)

Substituting the asymptotic forms for 212Pb and 220Rn in eq. (76) yields:

$\qquad\left(\lambda_{Pb}+\alpha_{Pb}-\lambda_{Ra}\right)\widetilde{n}_{Pb}=D_{Pb}\left(\frac{d^{2}\widetilde{n}_{Pb}}{dr^{2}}+\frac{1}{r}\frac{d\widetilde{n}_{Pb}}{dr}\right)+\lambda_{Rn}A_{Rn}K_{0}\left(\frac{r}{L_{Rn}}\right)$
(79)

where $A_{Rn}$ is given in (72). Using the definition of $L_{Pb}$ (10), this
can be written as:

$\left(\frac{d^{2}\widetilde{n}_{Pb}}{dr^{2}}+\frac{1}{r}\frac{d\widetilde{n}_{Pb}}{dr}\right)-\frac{1}{L_{Pb}^{2}}\widetilde{n}_{Pb}+\frac{\lambda_{Rn}}{D_{Pb}}A_{Rn}K_{0}\left(\frac{r}{L_{Rn}}\right)=0$
(80)

We attempt a solution of the form:

$\widetilde{n}_{Pb}\left(r\right)=A_{Pb}K_{0}\left(\frac{r}{L_{Rn}}\right)+B_{Pb}K_{0}\left(\frac{r}{L_{Pb}}\right)$
(81)

In order to proceed, we need to utilize some properties of the modified Bessel
functions. Again, we define $\xi=r/L$. Using the recursion relation (69), we
have:

$\displaystyle\frac{dK_{0}}{d\xi}=-K_{1}\left(\xi\right)$ (82)
$\displaystyle\frac{dK_{1}}{d\xi}=\frac{1}{\xi}K_{1}\left(\xi\right)-K_{2}\left(\xi\right)$
(83)

Another property of the modified Bessel functions (?) is:

$K_{n+1}\left(\xi\right)=K_{n-1}\left(\xi\right)+\frac{2n}{\xi}K_{n}\left(\xi\right)$
(84)

Thus:

$K_{2}\left(\xi\right)=K_{0}\left(\xi\right)+\frac{2}{\xi}K_{1}\left(\xi\right)$
(85)

Substituting $K_{2}\left(\xi\right)$ from (85) in (83) gives:

$\frac{dK_{1}}{d\xi}=-\left(K_{0}\left(\xi\right)+\frac{1}{\xi}K_{1}\left(\xi\right)\right)$
(86)

From (82) we have:

$\frac{d}{dr}K_{0}\left(\frac{r}{L}\right)=\frac{1}{L}\frac{d}{d\xi}K_{0}(\xi)=-\frac{1}{L}K_{1}\left(\frac{r}{L}\right)$
(87)

and from (86):

$\frac{d}{dr}K_{1}\left(\frac{r}{L}\right)=\frac{1}{L}\frac{d}{d\xi}K_{1}\left(\xi\right)=-\frac{1}{L}\left(K_{0}\left(\frac{r}{L}\right)+\frac{L}{r}K_{1}\left(\frac{r}{L}\right)\right)$
(88)

Using results (87) and (88), we get:

$\frac{d\widetilde{n}_{Pb}}{dr}=-\frac{A_{Pb}}{L_{Rn}}K_{1}\left(\frac{r}{L_{Rn}}\right)-\frac{B_{Pb}}{L_{Pb}}K_{1}\left(\frac{r}{L_{Pb}}\right)$
(89)

and:

$\displaystyle\quad\frac{d^{2}\widetilde{n}_{Pb}}{dr^{2}}=\frac{A_{Pb}}{L_{Rn}^{2}}\left(K_{0}\left(\frac{r}{L_{Rn}}\right)+\frac{L_{Rn}}{r}K_{1}\left(\frac{r}{L_{Rn}}\right)\right)+\frac{B_{Pb}}{L_{Pb}^{2}}\left(K_{0}\left(\frac{r}{L_{Pb}}\right)+\frac{L_{Pb}}{r}K_{1}\left(\frac{r}{L_{Pb}}\right)\right)$

resulting in:

$\frac{d^{2}\widetilde{n}_{Pb}}{dr^{2}}+\frac{1}{r}\frac{d\widetilde{n}_{Pb}}{dr}=\frac{A_{Pb}}{L_{Rn}^{2}}K_{0}\left(\frac{r}{L_{Rn}}\right)+\frac{B_{Pb}}{L_{Pb}^{2}}K_{0}\left(\frac{r}{L_{Pb}}\right)$
(91)

When (91) is inserted in eq. (80), the $K_{0}\left(r/L_{Pb}\right)$ terms
cancel out, leaving:

$\quad
A_{Pb}=\frac{L_{Rn}^{2}L_{Pb}^{2}}{L_{Rn}^{2}-L_{Pb}^{2}}\frac{\lambda_{Rn}}{D_{Pb}}A_{Rn}=\frac{L_{Rn}^{2}L_{Pb}^{2}}{L_{Rn}^{2}-L_{Pb}^{2}}\,\frac{\lambda_{Rn}}{D_{Pb}}\,\frac{P_{des}(Rn)\left(\Gamma_{Ra}^{src}(0)/l\right)}{2\pi
D_{Rn}\,(R_{0}/L_{Rn})K_{1}\left(R_{0}/L_{Rn}\right)}$ (92)

The coefficient $B_{Pb}$ is found from the boundary condition (77). Defining
$\widetilde{j}_{Pb}\left(r\right)=-D_{Pb}\frac{d}{dr}\widetilde{n}_{Pb}\left(r\right)$
and using (89), we have:

$\displaystyle\qquad 2\pi R_{0}\widetilde{j}_{Pb}(R_{0})$ $\displaystyle=2\pi
R_{0}D_{Pb}\left(\frac{A_{Pb}}{L_{Rn}}K_{1}\left(\frac{R_{0}}{L_{Rn}}\right)+\frac{B_{Pb}}{L_{Pb}}K_{1}\left(\frac{R_{0}}{L_{Pb}}\right)\right)$
(93)
$\displaystyle=\left(P_{des}^{eff}(Pb)-P_{des}(Rn)\right)\frac{\Gamma_{Ra}^{src}(0)}{l}$

yielding:

$\qquad
B_{Pb}=\frac{\left(P_{des}^{eff}(Pb)-P_{des}(Rn)\right)\left(\Gamma_{Ra}^{src}(0)/l\right)}{2\pi
D_{Pb}\,(R_{0}/L_{Pb})K_{1}(R_{0}/L_{Pb})}-A_{Pb}\frac{(R_{0}/L_{Rn})K_{1}(R_{0}/L_{Rn})}{(R_{0}/L_{Pb})K_{1}(R_{0}/L_{Pb})}$
(94)

For an ideal line source, $(R_{0}/L)K_{1}(R_{0}/L)\rightarrow 1$ giving:

$A_{Pb}^{line}=\frac{L_{Rn}^{2}L_{Pb}^{2}}{L_{Rn}^{2}-L_{Pb}^{2}}\frac{\lambda_{Rn}}{D_{Pb}}\frac{P_{des}(Rn)\left(\Gamma_{Ra}^{src}(0)/l\right)}{2\pi
D_{Rn}}$ (95)

$B_{Pb}^{line}=\frac{\left(P_{des}^{eff}(Pb)-P_{des}(Rn)\right)\left(\Gamma_{Ra}^{src}(0)/l\right)}{2\pi
D_{Pb}}-A_{Pb}$ (96)

Finally, we have:

$n_{Pb}^{asy}(r,t)=\left(A_{Pb}K_{0}\left(\frac{r}{L_{Rn}}\right)+B_{Pb}K_{0}\left(\frac{r}{L_{Pb}}\right)\right)e^{-\lambda_{Ra}t}$
(97)

###  A.3 212Bi

The diffusion-leakage equation for 212Bi for an infinite cylindrical source
is:

$\frac{\partial n_{Bi}}{\partial t}=\frac{D_{Bi}}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial n_{Bi}}{\partial
r}\right)-\lambda_{Bi}n_{Bi}-\alpha_{Bi}n_{Bi}+\lambda_{Pb}n_{Pb}$ (98)

Since 212Bi atoms are not emitted directly from the source, the boundary
condition is:

$\lim_{r\rightarrow R_{0}}2\pi rj_{Bi}(r,t)=0$ (99)

We look for an asymptotic solution of the form:
$n_{Bi}^{asy}(r,t)=\widetilde{n}_{Bi}(r)e^{-\lambda_{Ra}t}$, attempting:

$\widetilde{n}_{Bi}(r)=A_{Bi}K_{0}\left(\frac{r}{L_{Rn}}\right)+B_{Bi}K_{0}\left(\frac{r}{L_{Pb}}\right)+C_{Bi}K_{0}\left(\frac{r}{L_{Bi}}\right)$
(100)

Substituting the asymptotic form (100) in eq. (98) and in the boundary
condition (99), using the asymptotic form of the 212Pb solution (97) and the
expression for the effective 212Bi diffusion length (11), gives:

$\qquad
n_{Bi}^{asy}(r,t)=\left(A_{Bi}K_{0}\left(\frac{r}{L_{Rn}}\right)+B_{Bi}K_{0}\left(\frac{r}{L_{Pb}}\right)+C_{Bi}K_{0}\left(\frac{r}{L_{Bi}}\right)\right)e^{-\lambda_{Ra}t}$
(101)

For a cylindrical source of radius $R_{0}$:

$A_{Bi}=\frac{L_{Rn}^{2}L_{Bi}^{2}}{L_{Rn}^{2}-L_{Bi}^{2}}\frac{\lambda_{Pb}}{D_{Bi}}A_{Pb}$
(102)

$B_{Bi}=\frac{L_{Pb}^{2}L_{Bi}^{2}}{L_{Pb}^{2}-L_{Bi}^{2}}\frac{\lambda_{Pb}}{D_{Bi}}B_{Pb}$
(103)

$C_{Bi}=-\frac{(R_{0}/L_{Rn})K_{1}(R_{0}/L_{Rn})A_{Bi}+(R_{0}/L_{Pb})K_{1}(R_{0}/L_{Pb})B_{Bi}}{(R_{0}/L_{Bi})K_{1}(R_{0}/L_{Bi})}$
(104)

where $A_{Pb}$ and $B_{Pb}$ are given in (92) and (94). For an infinite line
source, $A_{Pb}$ and $B_{Pb}$ are replaced by their line source forms (95) and
(96) and:

$C_{Bi}^{line}=-(A_{Bi}+B_{Bi})$ (105)

##  Appendix B Discretization of the DL model equations in cylindrical
coordinates at the seed wall

The discretization of the DL equations for elements (cylindrical shells or
rings) touching the seed wall requires special attention, making use of the
appropriate boundary conditions. We demonstrate this in 1D, with practically
the same procedure for the 2D case. As stated in the main text, the radial
coordinate of the $i$-th cylindrical shell is given for the shell center:
$r_{i}=R_{0}+(i-\frac{1}{2})\Delta r$. For $i>1$ the diffusive term in the 1D
DL equations is discretized as follows:

$\displaystyle\qquad\frac{D}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial n}{\partial r}\right)_{r_{i}}\Longrightarrow$
$\displaystyle\qquad\frac{D}{r_{i}}\frac{\left(r\frac{\partial n}{\partial
r}\right)_{r_{i}+\Delta r/2}-\left(r\frac{\partial n}{\partial
r}\right)_{r_{i}-\Delta r/2}}{\Delta
r}=\frac{D}{r_{i}}\frac{(r_{i}+\frac{\Delta r}{2})\frac{n_{i+1}-n_{i}}{\Delta
r}-(r_{i}-\frac{\Delta r}{2})\frac{n_{i}-n_{i-1}}{\Delta r}}{\Delta r}=$
$\displaystyle\qquad\frac{D}{\Delta r^{2}}\left(\left(1+\frac{\Delta
r}{2r_{i}}\right)n_{i+1}+\left(1-\frac{\Delta
r}{2r_{i}}\right)n_{i-1}-2n_{i}\right)$ (106)

In the special case $i=1$, i.e., the shell immediately outside the source
surface with $r_{1}=R_{0}+\Delta r/2$, we use the boundary conditions at
$r\rightarrow R_{0}$. Since the current density is $j=-D\left(\partial
n/\partial r\right)$, we can write:

$\displaystyle\qquad\frac{D}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial n}{\partial r}\right)_{r_{1}}\Longrightarrow$
$\displaystyle\qquad\frac{D}{R_{0}+\Delta r/2}\frac{\left(R_{0}+\Delta
r\right)\left(\frac{\partial n}{\partial r}\right)_{R_{0}+\Delta
r}-R_{0}\left(\frac{\partial n}{\partial r}\right)_{R_{0}}}{\Delta
r}=\frac{D}{R_{0}+\Delta r/2}\frac{\left(R_{0}+\Delta
r\right)\frac{n_{2}-n_{1}}{\Delta r}+R_{0}\frac{j(R_{0},t)}{D}}{\Delta r}$

The source boundary conditions eq. (4)-(6) give:

$\displaystyle
R_{0}j_{Rn}(R_{0},t)=\frac{1}{2\pi}P_{des}(Rn)\frac{\Gamma_{Ra}^{src}(0)}{l}e^{-\lambda_{Ra}t}$
(108) $\displaystyle
R_{0}j_{Pb}(R_{0},t)=\frac{1}{2\pi}\left(P_{des}^{eff}(Pb)-P_{des}(Rn))\right)\frac{\Gamma_{Ra}^{src}(0)}{l}e^{-\lambda_{Ra}t}$
(109) $\displaystyle R_{0}j_{Bi}(R_{0},t)=0$ (110)

Therefore:

$\displaystyle\qquad\frac{D_{Rn}}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial n_{Rn}}{\partial r}\right)_{r_{1}}\Longrightarrow$
$\displaystyle\qquad\frac{D_{Rn}}{\Delta r^{2}}\left(\frac{1+\Delta
r/R_{0}}{1+\Delta
r/2R_{0}}\right)\left(n_{Rn,2}-n_{Rn,1}\right)+\frac{P_{des}(Rn)\left(\Gamma_{Ra}^{src}(0)/l\right)e^{-\lambda_{Ra}t}}{2\pi
R_{0}\Delta r\left(1+\Delta r/2R_{0}\right)}$

$\displaystyle\qquad\frac{D_{Pb}}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial n_{Pb}}{\partial r}\right)_{r_{1}}\Longrightarrow$
$\displaystyle\qquad\frac{D_{Pb}}{\Delta r^{2}}\left(\frac{1+\Delta
r/R_{0}}{1+\Delta
r/2R_{0}}\right)\left(n_{Pb,2}-n_{Pb,1}\right)+\frac{\left(P_{des}^{eff}(Pb)-P_{des}(Rn)\right)\left(\Gamma_{Ra}^{src}(0)/l\right)e^{-\lambda_{Ra}t}}{2\pi
R_{0}\Delta r\left(1+\Delta r/2R_{0}\right)}$

and

$\displaystyle\qquad\frac{D_{Bi}}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial n_{Bi}}{\partial
r}\right)_{r_{1}}\Longrightarrow\frac{D_{Bi}}{\Delta
r^{2}}\left(\frac{1+\Delta r/R_{0}}{1+\Delta
r/2R_{0}}\right)\left(n_{Bi,2}-n_{Bi,1}\right)$ (113)

## References

## References

  * [1]
  * [2] [] Abramowitz, M. & Stegun, I. A. (eds) (1972).  Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing, New York: Dover, chapter ”Modified Bessel Functions I and K” section 9.6, pp. 374–7.
  * [3]
  * [4] [] Arazi, L. (2020).  Diffusing alpha-emitters radiation therapy: approximate modeling of the macroscopic alpha particle dose of a point source, Physics in Medicine & Biology 65(1): 015015.
*#1
  * [5]
  * [6] [] Arazi, L., Cooks, T., Schmidt, M., Keisari, Y. & Kelson, I. (2007).  Treatment of solid tumors by interstitial release of recoiling short-lived alpha emitters., Phys Med Biol 52(16): 5025–42.
*#1
  * [7]
  * [8] [] Arazi, L., Cooks, T., Schmidt, M., Keisari, Y. & Kelson, I. (2010).  The treatment of solid tumors by alpha emitters released from 224 ra-loaded sources—internal dosimetry analysis, Phys Med Biol 55(4): 1203.
*#1
  * [9]
  * [10] [] Bellia, S., Feliciani, G., Duca, M., Monti, M., Turri, V., Sarnelli, A., Romeo, A., Kelson, I., Keisari, Y., Popovtzer, A., Ibrahim, T., Paganelli, G. & Stanganelli, I. (2019).  Clinical evidence of abscopal effect in cutaneous squamous cell carcinoma treated with diffusing alpha emitters radiation therapy: a case report, Journal of Contemporary Brachytherapy 11(5): 449–457.
*#1
  * [11]
  * [12] [] Confino, H., Hochman, I., Efrati, M., Schmidt, M., Umansky, V., Kelson, I. & Keisari, Y. (2015).  Tumor ablation by intratumoral ra-224-loaded wires induces anti-tumor immunity against experimental metastatic tumors, Cancer Immunology, Immunotherapy 64(2): 191–199.
*#1
  * [13]
  * [14] [] Confino, H., Schmidt, M., Efrati, M., Hochman, I., Umansky, V., Kelson, I. & Keisari, Y. (2016).  Inhibition of mouse breast adenocarcinoma growth by ablation with intratumoral alpha-irradiation combined with inhibitors of immunosuppression and cpg, Cancer Immunology, Immunotherapy 65(10): 1149–1158.
*#1
  * [15]
  * [16] [] Cooks, T., Arazi, L., Efrati, M., Schmidt, M., Marshak, G., Kelson, I. & Keisari, Y. (2009b).  Interstitial wires releasing diffusing alpha emitters combined with chemotherapy improved local tumor control and survival in squamous cell carcinoma-bearing mice, Cancer 115(8): 1791–1801.
*#1
  * [17]
  * [18] [] Cooks, T., Arazi, L., Schmidt, M., Marshak, G., Kelson, I. & Keisari, Y. (2008).  Growth retardation and destruction of experimental squamous cell carcinoma by interstitial radioactive wires releasing diffusing alpha-emitting atoms., Int J Cancer 122(7): 1657–64.
*#1
  * [19]
  * [20] [] Cooks, T., Schmidt, M., Bittan, H., Lazarov, E., Arazi, L., Kelson, I. & Keisari, Y. (2009a).  Local control of lung derived tumors by diffusing alpha-emitting atoms released from intratumoral wires loaded with radium-224, International Journal of Radiation Oncology * Biology * Physics 74(3): 966–973.
*#1
  * [21]
  * [22] [] Cooks, T., Tal, M., Raab, S., Efrati, M., Reitkopf, S., Lazarov, E., Etzyoni, R., Schmidt, M., Arazi, L., Kelson, I. & Keisari, Y. (2012).  Intratumoral 224ra-loaded wires spread alpha-emitters inside solid human tumors in athymic mice achieving tumor control, Anticancer Research 32(12): 5315–5321.
*#1
  * [23]
  * [24] [] Domankevich, V., Cohen, A., Efrati, M., Schmidt, M., Rammensee, H.-G., Nair, S. S., Tewari, A., Kelson, I. & Keisari, Y. (2019).  Combining alpha radiation-based brachytherapy with immunomodulators promotes complete tumor regression in mice via tumor-specific long-term immune response, Cancer Immunology, Immunotherapy 68(12): 1949–1958.
*#1
  * [25]
  * [26] [] Domankevich, V., Efrati, M., Schmidt, M., Glikson, E., Mansour, F., Shai, A., Cohen, A., Zilberstein, Y., Flaisher, E., Galalae, R., Kelson, I. & Keisari, Y. (2020).  Rig-1-like receptor activation synergizes with intratumoral alpha radiation to induce pancreatic tumor rejection, triple-negative breast metastases clearance, and antitumor immune memory in mice, Frontiers in Oncology 10: 990.
*#1
  * [27]
  * [28] [] Horev-Drori, G., Cooks, T., Bittan, H., Lazarov, E., Schmidt, M., Arazi, L., Efrati, M., Kelson, I. & Keisari, Y. (2012).  Local control of experimental malignant pancreatic tumors by treatment with a combination of chemotherapy and intratumoral 224radium-loaded wires releasing alpha-emitting atoms, Translational Research 159(1): 32–41.
*#1
  * [29]
  * [30] [] Keisari, Y., Hochman, I., Confino, H., Korenstein, R. & Kelson, I. (2014).  Activation of local and systemic anti-tumor immune responses by ablation of solid tumors with intratumoral electrochemical or alpha radiation treatments, Cancer Immunology, Immunotherapy 63(1): 1–9.
*#1
  * [31]
  * [32] [] Keisari, Y. & Kelson, I. (2021).  The potentiation of anti-tumor immunity by tumor abolition with alpha particles, protons, or carbon ion radiation and its enforcement by combination with immunoadjuvants or inhibitors of immune suppressor cells and checkpoint molecules, Cells 10(2).
*#1
  * [33]
  * [34] [] Keisari, Y., Popovtzer, A. & Kelson, I. (2020).  Effective treatment of metastatic cancer by an innovative intratumoral alpha particle-mediated radiotherapy in combination with immunotherapy: A short review, Journal of Physics: Conference Series 1662: 012016.
*#1
  * [35]
  * [36] [] Kis, T. (2021).  Tridiagonal matrix algorithm (thomas alg.) (tridiagonal).  Accessed: September 28, 2021.
*#1
  * [37]
  * [38] [] Lazarov, E., Arazi, L., Efrati, M., Cooks, T., Schmidt, M., Keisari, Y. & Kelson, I. (2011).  Comparative in vitro microdosimetric study of murine- and human-derived cancer cells exposed to alpha particles, Radiation Research 177(3): 280–287.
*#1
  * [39]
  * [40] [] Milrot, E., Jackman, A., Flescher, E., Gonen, P., Kelson, I., Keisari, Y. & Sherman, L. (2013).  Enhanced killing of cervical cancer cells by combinations of methyl jasmonate with cisplatin, x or alpha radiation, Investigational New Drugs 31(2): 333–344.
*#1
  * [41]
  * [42] [] Popovtzer, A., Rosenfeld, E., Mizrachy, A., Bellia, R. S., Ben-Hur, R., Feliciani, G., Sarnelli, A., Arazi, L., Deutsch, L., Kelson, I. & Keisari, Y. (2019).  Initial safety and tumor control results from a “first-in-human” prospective trial evaluating a novel alpha-emitting radionuclide for the treatment of locally advanced recurrent squamous cell carcinomas of the skin and head and neck, Accepted for publication in the International Journal of Radiation Oncology Biology Physics .
  * [43]
  * [44] [] Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (2007).  Numerical recipes 3rd edition: The art of scientific computing, Cambridge university press.
  * [45]
  * [46] [] Reitkopf-Brodutch, S., Confino, H., Schmidt, M., Cooks, T., Efrati, M., Arazi, L., Rath-Wolfson, L., Marshak, G., Kelson, I. & Keisari, Y. (2015).  Ablation of experimental colon cancer by intratumoral 224radium-loaded wires is mediated by alpha particles released from atoms which spread in the tumor and can be augmented by chemotherapy, International Journal of Radiation Biology 91(2): 179–186.
  * [47]