Comment on “Reproducibility study of Monte Carlo simulations for nanoparticle dose enhancement and biological modeling of cell survival curves” by Velten et al.

This comment highlights two methodological issues with the recent article by Velten et al. [ Biomed Phys Eng Express 2023;9:045004]


Introduction
In their recent article, Velten et al. [1] replicated the physical simulation and biological modeling of previously published research [2] on 50 nm gold nanoparticles (GNPs) irradiated with different photon spectra or with protons in a spread-out Bragg peak (SOBP).For the biological modeling, they used a variant of the local effect model (LEM), which they called the "LEM approach in 2D" and which is based on the assumption of an infinitesimally thin cell.This comment is to point out that (a) this "LEM approach in 2D" leads to biased results, and (b) there is a conceptional problem with the quantity referred to as a "macroscopic dose".

Concerns regarding the "LEM approach in 2D"
The LEM in its original form [3] is based on the assumption that the survival rate S of a cell line at an absorbed dose D administered by low linear energy transfer (LET) irradiation is given by () =  −〈  ()〉 (1) where 〈  ()〉 is the expected number of lethal lesions in the nucleus (at this dose and for this low-LET radiation quality).For low-LET radiation, the statistical expectation of the energy imparted in a small subvolume of the nucleus is assumed to be uniformly distributed across the nucleus.
Similarly, the probability that a lesion will occur in such a small subvolume is also assumed to be uniformly distributed over the nucleus.Therefore, if the irradiation is such that the expected value of the energy imparted varies across the dimensions of the cell nucleus, as for ion beams, the expected value of the number of lethal lesions is given by Eq. (2).

〈𝑁(𝐷; 𝛼, 𝛽)〉 = ∭
((, , ); , ) Here,   is the volume of the nucleus and ,  are the parameters of the linear-quadratic (LQ) model describing cell survival as a function of dose for uniform irradiation with low-LET radiation.〈(; , )〉 is thus the integral of the spatial density of the expected lesions.
Although not explicitly stated, the description in the Materials and Methods section suggests that the further steps of their procedure were as follows.First, a number of points equal to the expected number of nanoparticles undergoing a photon or proton interaction were randomly selected in the infinitesimally thin cell.Then the dose kernels obtained from their simulations were placed at these positions to obtain the two-dimensional dose distribution, which was then used in Eq. ( 3) to obtain the expected number of lethal lesions.
If this interpretation of the procedure is correct and if the dose kernels shown in Supplementary Fig. 5 of Velten et al. (2023) were used, this leads to an overestimation of the number of lesions and an underestimation of cell survival.This is because using a cross-section of the 3D dose profile is equivalent to assuming that the cell has the shape of a cylinder in which the dose profile is constant along the direction of the cylinder axis.(As is approximately true for ion beams with energies much higher than the Bragg peak energy.)But then, instead of a nanoparticle undergoing an interaction, which is approximately a point source, one would have a line source.This would result in a much higher average dose in the cell nucleus than would be the case for nanoparticles.
Using a cut through the 3D dose map instead of an Abel transform from 3D to 2D means that the dose profile has both wrong magnitude and wrong dependence on radial distance from the GNP.Transforming the 3D dose distribution to 2D is equivalent to averaging over the third dimension.With a 1/r² dose dependence in 3D, the transformed distribution has a 1/r dependence, and the peak dose at the surface of the GNP is reduced by a factor proportional to the ratio of GNP diameter to the thickness of the nucleus.The proportionality constant is of the order of .For a GNP radius of 25 nm and a nucleus diameter of 8 µm, the reduction factor is thus about 0.01, which means that the dose at the surface of the GNP is 3 Gy for the 83 keV photons (which produce the highest surface dose, cf.Table S2 in [1]) instead of 300 Gy.
Since nothing of the preceding is mentioned in the work of Velten et al., it must be assumed that a central cut through the 3D dose distribution was used and that the results are overestimated as explained in the following.

An illustrative example: a single nanoparticle
To illustrate the argument above, consider the case where the GNPs are inside the nucleus and only one GNP undergoes an interaction.Assume further that the extra dose due to an ionization in the GNP varies with the inverse square of the radial distance from the GNP center.Then the dose inside this cell nucleus is essentially equal to D0 everywhere except inside a sphere of radius rm around the GNP, where the dose significantly exceeds D0 and is given by Eq. ( 4).
Here  1 is the excess dose at the surface of the spherical GNP with radius   , and the origin of the coordinate system is at the center of the GNP.Assume further that the GNP is located such that aforementioned sphere of radius rm is completely inside the cell nucleus.
For simplicity, we further assume that the dose at the surface of the GNP is below the threshold dose for the transition from the LQ dose dependence of cell survival to linear high-dose dependence.Then the expected number of lesions for the 3D case is obtained as For the 2D case, assume that the nucleus has a cross section An and that the GNP is located such that the circle of radius rm around the GNP is completely inside the nucleus.
Then the dose in the nucleus is D0 outside the circle, while the dose inside the circle is given by Eq. ( 4), where the origin of the coordinate system is again at the center of the GNP.The expected number of lesions is then obtained as The ratio  3 between the third summand in Eq. ( 6) and the third summand in Eq. ( 5) is where hn = Vn /An, that is, the average thickness of the nucleus along the third dimension.
The corresponding ratio  4 between the fourth summands in Eq. ( 6) and Eq. ( 5) is In Fig. 1 of Velten et al., the dose around a GNP undergoing an interaction falls to below 0.5 Gy within 250 nm of the GNP center.For the photon irradiation, the decrease between the surface of the GNP and this distance is even more than two orders of magnitude, since the gold Mshell Auger electrons are completely stopped in this range.Thus, 250 nm appears to be a reasonable lower limit for   .If the sphere and the circle both have a diameter of 8 µm, as considered by Velten et al. [1], then hn is about 5.3 µm.Using these values in Eqs.(7) and (8),  3 and  4 are obtained as approximately 30 and 60, respectively.This means that the number of lesions originating from the GNP is overestimated by orders of magnitude in this case when the "LEM approach in 2D" is used instead of the conventional LEM.

Adaptation of the example to many nanoparticles
The case of many GNPs instead of just one is more complex.Since the GNPs are assumed to be uniformly distributed in the cell nucleus, a modified version of Eq. ( 4) can be used for the dose in the vicinity of the GNP: Here   is the background dose, which also contains the contributions from the other GNPs undergoing an ionizing interaction.As is shown in Fig. S1.1 (Supplement 1), this background dose varies over the cell nucleus by up to a factor of 2 (in the 3D case).The mean values for the 3D and 2D cases (cf.Supplement 1) are given by Eqs. ( 10) and (11).
Here,  0 is the dose in the absence of GNPs, Rn is the radius of the nucleus, and   =  0 is the number of GNPs undergoing an ionizing interaction.N is the total number of GNPs in the nucleus, and p is the probability per dose of a GNP undergoing an ionizing interaction.
Neglecting the variation of   and using only the mean value  ̅  , the necessary adjustments to Eqs. ( 4) and ( 5) are to replace  0 by  ̅  ,  1 by    1 , and  1 ² by    1 ², which leads to Eqs. ( 12) and (13).
Here it was taken into account, that the volumes or areas of the GNPs must be excluded when calculating the integrals in Eqs. ( 2) and (3).The corrected values of nucleus volume   ′ and area   ′ are given in Eq. ( 14).
It is important to note that also the contributions from the background dose are different between the 3D and 2D cases.Generally,  ̅ ,2 is larger than  ̅ ,3 if the number of NPs undergoing an interaction is higher than 1 because then  2 is smaller than  3 (cf.Eq. (S1.16) in Supplement 1).

Numerical example
It may be useful to consider concrete values by using the parameters from the work of Velten et al., namely  = 0.019 Gy -1 ,  = 0.052 Gy -2 , rs = 25 nm, and Rn = 4 µm.
Assume that 1000 GNPs undergo an ionizing interaction.(This corresponds approximately to the case of 80 keV photons at a dose without GNPs of D0 = 1, cf.Table 1 in [1]).With an assumed surface dose of the GNP of D1 = 100 Gy,  ̅ ,3 and  ̅ ,2 are about 8.6 Gy and 26.5 Gy, respectively.The resulting predicted number of lesions is 5.30 and 61.3, respectively.Here, it was taken into account that the background dose in the 2D case exceeds the threshold dose of 20 Gy for the transition from the linear-quadratic to a linear dose dependence.(See Eq. (S2.9) in Supplement 2 for the expression replacing Eq. ( 13) in this case.)That is, the 2D approach predicts more than 10 times more lesions than the 3D approach in this case.
The contribution of the background dose to the predicted lesions is 4.02 in the 3D LEM, while 1.17 and 0.11 lesions are coming from the third and fourth terms in Eq. ( 12), respectively.In the LEM in 2D case, the background dose accounts for 34.7 lesions, and the dose peaks around GNPs contribute 26.6 lesions.These numbers corroborate that a significant overestimation of the number of lesions occurs in the "LEM in 2D" when the 3D dose kernels are used.Of course, the concrete numbers are only valid for the specific case and the assumptions made in the derivation, such as the 1/r² dose dependence.
The number of interacting GNPs increases proportional to the incident fluence (or dose without GNPs).At higher numbers of interacting GNPs, the contribution of the background dose grows faster than the LEM-specific terms in Eqs. ( 12) and (13) relating to dose non-uniformity (cf.Supplement 3).However, the overall overestimation of the number of lesions remains, as can be seen in Fig. 1.
The red curve in Fig. 1 relates to the LEM-specific third term of Eqs. ( 12) and (13).Here, the ratio between the 2D version and the 3D version of the LEM is over a factor of 10 over most of the range of dose without GNP.(The ratio of the fourth terms rapidly drops to zero when the background dose of the 2D case exceeds the threshold dose for transition from the linear-quadratic to a linear dose dependence.)The ratio of the total number of lesions from the 2D and 3D models settles around 5 in Fig. 1.
Fig. 1: Ratios R3 and R4 of the third and fourth terms in Eq. ( 13) to the third and fourth terms in Eq. ( 12) as a function of dose without GNPs (red and blue curve) for the model parameters considered in the text.The black curve shows the ratio of the total number of lesions predicted by Eq. ( 13) to the value obtained with Eq. ( 12).
How this overestimation of the number of lesions affects the predicted survival curves can be seen in Fig. 2, where the upper panel refers to the case considered so far, namely that 1000 GNPs undergo an ionizing interaction at a dose of 1 Gy.The bottom panel refers to the case of only 100 ionizations in GNPs at 1 Gy.This could occur due to different energies of the incident radiation or when lower concentrations of GNPs are considered.In both cases, it can be seen in Fig. 2 that the LEM in 2D predicts a much lower surviving fraction than the 3D version of the LEM.

Concerns regarding "macroscopic dose" estimates
Absorbed dose is a point quantity.Therefore, it is reasonable to assume that what Velten et al. and Lin et al. [2] refer to as the "macroscopic dose" is the macroscopic average dose.More precisely, it appears to be the macroscopic average dose to water in the absence of GNPs.
It is not clear from either work whether an isolated cell loaded with GNPs and surrounded by cells without GNPs is being considered or a tumor loaded with GNPs.In the latter case, all but the cells at the tumor's surface would be in an environment with the same gold concentration as the cell under consideration itself.In a more realistic scenario, where the tissue adjacent to the tumor also contains GNPs, the preceding statement applies to all tumor cells.
The macroscopic average dose in a piece of tissue with a given concentration of gold in the form of GNPs is enhanced by the increased absorption of gold compared to the elements that tissue is composed of.For the photon spectra shown by Velten et al., the macroscopic average dose can be estimated from the mass concentration of gold and from literature data on mass energy-absorption coefficients [4].
For the photon irradiations, 1.410 5 GNPs of 50 nm diameter were considered, which have a total volume of about 9 µm³.This corresponds to 3.4 % of the volume of the nucleus (8 µm diameter sphere) and 0.7 % of the volume of the cell (13.5 µm diameter).The corresponding mass fractions of gold are 660 g/kg and 137 g/kg.
The gold-to-water ratio of the mass energy-absorption coefficients for the monoenergetic photon energies ranges from 11 (250 keV) to 145 (50 keV).If the weights of the differrent photon energies in the mock 250 kVp spectra are used (given in Table S1 of [1]), this ratio is between 58 and 87.Therefore, the macroscopic average dose under secondary particle equilibrium for the case of the GNPs uniformly distributed over the cells is expected to be higher.This dose is expected to be higher by a factor between 2.4 (250 keV) and 20.7 (50 keV)or between 8.8 and 12.8 for the mock 250 kVp spectrathan for the case without GNPs.
These factors are even higher than the radiation enhancement factors (ratio between the doses leading to the same survival rate) seen in the left part of Fig. 5 in [1].The reason for this is as follows.In the approach of Velten et al. the proportion of the energy of electrons produced in GNPs, which is imparted outside the cell, is not compensated by energy imparted in a realistic irradiation scenario by electrons from radiation interactions in GNPs located in neighboring cells.
In the other two cases, where the GNPs are located only in the cell nucleus and only in the cytoplasm, the dose inside the nucleus differs from that in the cytoplasm.In fact, there is a dose gradient at the interface between the two regions.In the simple geometry of concentric spheres, the variation of the absorbed dose with radial distance from the center of the spheres could be calculated from the data shown in Fig. 1 of [1].
In any case, the "macroscopic dose" plotted on the x-axes in Fig. 5 and Fig. S6 of [1] is not the macroscopically averaged dose for the curves with GNPs and should be called differently.For instance, the incident photon fluence could be used instead.Or it should be clearly stated that it is the dose in the absence of GNPs as has been done with Fig. 1 and Fig. 2 of this paper.(For comparing the outcomes of the two models, the use of dose without GNPs as independent variable was a necessity.)

Conclusions
It is understood that Velten et al. [1] intended to reproduce the work of Lin et al. [2] and add some new insights into the contributions of photons slightly above the gold K-edge.The author fully endorses their recommendations to report more details on the photon spectra used in publications.Such details are a prerequisite for any attempt to check whether other codes give the same results.
Insufficient details on the methods used make it difficult for other researchers to assess or verify the results.Unfortunately, like Lin et al. before them, Velten et al. also fail to provide sufficient details on their modifications to the LEM or to justify why this differing approach should give correct predictions.From the analysis given in this Comment, it appears that the "LEM approach in 2D" as used by Velten et al. can lead to survival rates that are underestimated by orders of magnitude.
Another recommendation made here is to only use welldefined terms or to clearly state what is meant, for instance, by "macroscopic dose".Plotting results for cell survival in the presence of GNPs as a function of the macroscopic average dose in their absence is justified for low gold concentrations.For gold concentrations as high as 660 g/kg and 137 g/kg (which appear unrealistic compared to values reported in radiobiological studies [5]), this approach is misleading.The use of the average dose with GNPs seems better suited to disentangle the effects of dose enhancement and "radio-enhancement" [6], that is, the increase in average dose versus local dose (that affects cell survival predicted by the LEM).

Supplement 1
Hans Rabus Physikalisch-Technische Bundesanstalt, Berlin, Germany Email: hans.rabus@ptb.deIn this Supplement, the mathematical expressions for the average dose in the cell nucleus and the background dose in the vicinity of a nanoparticle (NP) are derived for the case that • the NPs are uniformly distributed inside the nucleus, • the extra dose from an interaction in the NP varies with the inverse square of the radial distance from the center of the NP.
The nucleus is assumed to be a sphere and a circle in the 3D model and 2D model, respectively.Both sphere and circle have a radius of Rn.Then the number densities of NPs undergoing an ionizing interaction in the two cases,  3 and  2 , are given by Eq. (S1.1).
Here, N is the total number of NPs in the nucleus,  0 is the dose in the absence of NPs, and p is the probability of a NP undergoing a photon interaction for a photon fluence producing a dose of 1 Gy in the absence of NPs.(p corresponds to the parameter Ptot in the work of Velten et al. [1].) The additional dose   from a single NP undergoing an ionizing interaction is given by Eq. (S1.2).
Here  1 is the dose contribution at the surface of the NP due to an ionization in the NP, and   is the radius of the spherical NP.The average doses at a point in the nucleus are given by Eqs.(S1.3) and (S1.4) for the 3D case and 2D case, respectively, where  0 is the dose in the absence of NPs.
The support of the inner integrals appearing on the right-hand sides of Eqs.(S1.3) and and (S1.4) is between   and the length  of the chord from the point where the dose is evaluated to the surface of the nucleus along the direction specified by the polar angle θ and plane angle φ in the 3D case and 2D case, respectively.
The functional form of the expression for the length of the chord is given by Eq. (S1.7) (Fig. S1.1; in the 2D case,  must be replaced by ).
Here, r is radial distance of the point at which the dose is evaluated from the center of the nucleus.In the 3D case, the integral can be calculated analytically, and the last factor in Eq. (S1.5) is obtained as the expression given in Eq. (S1.8).
The mean value of  3 is given by Where in the last integral     ⁄ has been substituted with x.This integral can be calculated analytically and gives a value of 1/3 (see Appendix, page S1-5), so that finally the mean dose in the 3D case is obtained by Eq. (S1.10).

𝐷 ̅ 𝑎
The corresponding factor for the 2D case cannot be obtained analytically.The respective data shown as red curve in Fig. 1 were determined by Monte Carlo integration using an ad-hoc GNU data language (GDL) script assuming the parameter values rs = 25 nm and Rn = 4 µm.This script was also used for calculating the contribution to  ̅ 2 which is independent of rs.The resulting expressions for  ̅ 2 and  ̅ ,2 are given in Eqs.(S1.12) and (S1.13).The background dose in the vicinity of an NP undergoing an ionizing interaction produced by the other interacting NPs can be obtained in a similar way by considering only NP locations at radial distances exceeding those of the radius of a sphere of volume   /  (or of a circle of area   /  ) centered at the considered NP. Here is the number of NPs undergoing an ionizing interaction, and the radius of the sphere,  3 , and the radius of the circle,  2 , are given by:  Appendix: Evaluation of the last integral in Eq. (S1.9) The integral can be rewritten as follows: ∫ (  −  3 ) ln ( 1 +  1 −  )

Hans Rabus
Physikalisch-Technische Bundesanstalt, Berlin, Germany Email: hans.rabus@ptb.de In this Supplement, the mathematical expressions for the modifications of Eqs. ( 12) and (13) in the main document are derived for the case that the dose dependence of the number of induced lesions changes from linearquadratic to a linear when the dose D exceeds a threshold Dt.The derivation refers to spherical nanoparticles (NPs) located inside a spherical nucleus.
The linear dose dependence of the induced lesions at high doses is given by [1]: where α and β are the parameters of the linear-quadratic dose dependence at doses below Dt.Eq. (S2.1) can be conveniently rewritten as follows: In the vicinity of an NP, the dose is assumed to have the following dependence on radial distance r from the center of the NP: where Db is the background dose due to interactions in water and the other NPs undergoing an ionizing interaction, which was derived in Supplement 1. D1 is the extra dose at the NP surface due to an ionizing interaction in the NP, and ts is the radius of the NP.

Case 1: Db  Dt
In case that the background dose Db exceeds the threshold dose Dt, the dose around an NP is always larger than Dt as well.For a dose dependence as given in Eq. (S2.3) the following expression is obtained for the predicted number of lesions in the 3D case: Here Ni is the number of NPs undergoing an ionizing interaction, Vn is the volume of the nucleus, Vs is the volume of a NP, and N is the total number of NPs in the nucleus.The upper limit of the integral, r3, is the radius of a sphere having a volume Vn/Ni.
The corresponding expression for the 2D case is where An is the area of the 2D-nucleus, As is the cross-sectional area of a NP, and the upper limit of the integral, r2, is the radius of a circle of area An/Ni.

S2-2
Evaluating the integrals in Eqs.(S2.4) and (S2.If this threshold value is larger than r2 or r3 in the 2D or 3D model, respectively, then the dose inside the corresponding sphere or circle around an interacting NP always exceeds the threshold dose.Therefore, Eqs.(S2.8) are (S2.9)also describe this case.However, to distinguish between the lesions due to the background dose alone and the contributions due to dose-non-uniformity, it is convenient to rewrite the expressions as in Eqs.(S2.11) and (S2.12).

Fig. 2 :
Fig. 2: Survival curves obtained with the 3D and 2D versions of the LEM for the model parameters used by Velten et al. [1] and for (a) 1000 and (b) 100 GNPs undergoing an ionization at a dose without GNPs of 1 Gy.

Fig. S1. 1 :
Fig. S1.1: Illustration of the determination of the mean chord from a point to the surface of a sphere (or to the perimeter of a circle in the 2D case).
simply verified by calculating the derivative, the primitives of   ln  are ∫   ln   = Comment on "Reproducibility study of Monte Carlo simulations for nanoparticle dose enhancement and biological modeling of cell survival curves" by Velten et al.
6)gives Eqs.(S2.6) and (S2.7), respectively.Using the definitions of r2 and r3, Eqs.(S2.8) and (S2.9) are obtained as the final expressions.If the background dose Db is below the threshold dose Dt, then the dose in the nucleus exceeds Dt only within radial distances from a NP of below a threshold value rt given by Eq. (S2.10).