TG-43 comparison with egs_brachy simulations showing discrepancies near phantom edge

[vahx129]

Describe the bug I am not 100% confident this is a bug maybe so much as me trying to perform this validation incorrectly. But essentially, I'm trying to validate egs_brachy dose distributions in voxelized water phantom geometries by comparing them to TG-43 calculated doses from MIM which we use as our TPS. When I do the comparisons for post implant CTs from actual patient treatments, the ones where the seeds are relatively centrally placed in the axial plane of the CT image agree well with egs_brachy simulations. In the cases where the seeds are closer to the phantom edges, I'm seeing a lot more disagreement which I don't expect. It makes me think that maybe TG-43 scatter conditions aren't being faithfully reproduced even though I'm embedding my phantom geometry into a much larger box of water to approximate the infinite water conditions of TG-43.

To Reproduce A precise description of the steps to reproduce the behavior:

1. I convert a CT scan of a patient into an .egsphant geometry with a ramp that converts all densities to water and resize the voxels to be 0.5 mm in all directions.
2. I use a MATLAB script to place the Cs-131 seeds in the same positions as they were localized in the patient's post-implant scan and create an egs_brachy input file out of it that calls the .egsphant file (which is all water), embeds the seeds into their accurate treatment locations in that .egsphant file, and then embed the "seeds in egsphant" geometry into a very large box of water.
3. Run the simulation to create .3ddose files from egs_brachy and use a normalization to convert units to Gy based on mean lifetime and a TG-43 simulation of a single Cs-131 seed (cross-checked against the results from Carleton and came within about 0.5% of their normalization number). The seed strength in U is accounted for here as well since we only use uniform strength seeds when treating H&N patients.
4. Compare the 3ddose distribution in units of Gy to the DICOM RT-Dose file generated by the MIM TPS after seeds were localized in the post-implant CT scan.
5. Run various Gamma analyses to see that regions closer to the .egsphant edge have worse agreement with the TG-43 calculated doses than elsewhere.

Expected behavior I expect that my TG-43 TPS calculated dose distributions and Monte Carlo simulations, which are performed in a pure voxelized water phantom geometry itself embedded in a giant surrounding box of water to simulate infinite water conditions, should agree regardless of how close the seeds are to the edge of the voxelized portion of the phantom geometry. I can't say for certain that this isn't intended behavior based on how I'm running these simulations but when I look at the tracks in the "final" geometry of my simulations they look reasonable but since I'm only exporting the .3ddose of just the voxelized water phantom geometry to compare to my TPS dose, I'm wondering why the higher levels of disagreement are showing up so starkly near the phantom edge and not elsewhere. In cases where the seeds are centrally located within the phantom geometry, this doesn't seem to pop up at all.

Screenshots TPS Dose for one Axial Slice of clinical seed arrangement based on TG-43

egs_brachy dose for same axial slice with same clinical seed arrangement in a full water phantom

Global Gamma 2%/2mm dose comparison between TG-43 based calculation and egs_brachy in water phantom geometry calculation illustrating a lot of disagreement towards the edge of the phantom

egs_brachy tracks when visualizing just voxelized water phantom + seed geometry (same perspective as above images)

egs_brachy tracks when visualizing entire geometry which is seeds embedded in a voxelized water phantom matching resolution of CT scan itself embedded in a larger box of water

TPS Dose for a more centrally located set of seeds

egs_brachy dose for a more centrally located set of seeds

More favorable Gamma comparison between the previous two dose distributions

Operating system

• [X] Linux
• [ ] Windows
• [ ] macOS
• [ ] Other:

EGSnrc version egs_brachy 2024

Additional context I can share .egsphant, ramp, and input files if needed but am also still half of the mindset that I might just not be running these simulations properly. Either way, help would be hugely appreciated!

*Edited because I accidentally pasted the same screenshot twice of egs_brachy dose distribution instead of actual TPS one

[mchamberland]

That's an interesting one! Nothing immediately jumps to me.

Please share everything, so I can take a closer look.

If you need my email:

marc.chamberland@uvmhealth.org

[mchamberland]

Also, instead of gamma, could you look at the distribution of the dose differences in units of the statistical uncertainties in the egs_brachy voxels?

In other words, calculate (egs_brachy dose - TPS dose) / (absolute statistical uncertainty) on a voxel-per-voxel basis. Do it separately for both the central and edge case. You'd expect the two distributions to be roughly normally distributed and centred at 0.

[vahx129]

Thanks for the prompt response Marc! Hope you're doing well. For starters, here is a link to my relevant EGSnrc files which shouldn't contain PHI. I will work on anonymizing the two cases if requested for CT and DICOM-RT dose date if needed and also RT plan data to verify seed positioning (I'm pretty sure that isn't the issue here though).

Link to egs_brachy files: https://1drv.ms/f/s!ArqLIYhAstNSvuMlDY8JydPkzB7WDA?e=SxFr4m

At a glance of the relative uncertainties, I'm not seeing any indication that there is additional statistical uncertainty corresponding to the locations where the dose comparison appears to be worse. See the following:

Slice of dose and relative error side by side.

Corresponding Gamma Map

I should also mention that I do crop the CT in the slice or superior/inferior direction in order to save on simulation time which doesn't show up in the images I provided but we still have about 5-8 cm of thickness in that direction (and ostensibly much more water since it's embedded in the water box).

[mchamberland]

Thanks, I'll try to take a closer look later this week.

I still would like to see the distribution (histogram) of (egs_brachy dose - TPS dose) / (absolute uncertainty) at the voxel-by-voxel level. Roughly 95% should be within +/- 2. If that is the case, then the gamma analysis itself could be the issue.

[vahx129]

Still working on the voxel-by-voxel comparison. Because the egs_brachy 3ddose and TPS DICOM-RT Dose matrices are different sizes, I am going to try to write something up to crop down the TPS dose distribution to an identically sized matrix to do the subtraction appropriately. I won't go too crazy with interpolation but bear with me a bit longer as I try to put it together. It's been a while since I've worked with 3D dose distributions in MATLAB so I'm a little rusty. Thanks!

[mchamberland]

Two things:

• The egsinp you shared has the big box made of air, not water. If this is the correct egsinp, then there's your problem.

• You're not really doing a fair TG-43 comparison because you're using the normal run mode. You should be using superposition mode so that only one seed is active at a time, to properly simulate TG-43 conditions. The normal run mode includes interseed effects.

[vahx129]

I may have have sent over an incorrect input where I was messing with the "world box" materials (it is usually set to water). The normal run mode is a great catch though because that was an underlying issue with my egs_brachy input generation script using an incorrect template file.

When I try to manually convert the input to a superposition run, I do the following:

• Add "type = EGS_ASwitchedEnvelope" to my phantom_and_seeds geometry
• Change run type to "superposition"
• Add "source envelope geometry = phantom_and_seeds" near bottom of geometry definition

When I do this, I'm getting the following error:

Thinking this is hopefully something silly I'm messing up or forgetting on my end here but any thoughts as to what might be causing that error? I've made the input file called "HN6_AWat_5mmXYZ_2cm_Crop_Artifact_1e9_WB_SP.egsinp" available in the root directory of the folder I shared earlier.

Thanks a lot!

[vahx129]

Never mind it turned out Notepad++ was doing something strange with newline formatting. Have had this issue before and just switched to Notepad to rectify. Sorry to spam! Will let you know results of rerun simulation using superposition mode. Thanks!

[ftessier]

In Notepad++ you can set the end-of-line characters(s) under Edit $\rightarrow$ EOL conversion. Or, you can use VS Code 😉

[vahx129]

OK so with the new simulation run in superposition mode with a "world_water" box, I'm still seeing similar results. For the histogram analysis, I'm not sure if these units square with what Marc requested earlier but I essentially crop and interpolate my TPS dose data to the same size and resolution as the .3ddose matrix (1 mm voxels to 0.5 mm voxels) and then directly subtract the two matrices before normalizing by relative error multiplied by voxel dose within each voxel.

The resulting histogram of just (egs_brachy dose - TPS dose) looks like the following:

The resulting histogram of (egs_brachy dose - TPS dose) / (relative uncertainty * dose in each voxel) looks like the following:

Here are the gamma analyses compared for my old input (water box world but run in normal mode) and the new input (water box world but run in superposition)

Old Result:

New Result:

[mchamberland]

In this situation, I would take the TPS dose as having no uncertainty, so the denominator should just be (egs_brachy voxel relative uncertainty * egs_brachy voxel dose). Basically, you're trying to see if the dose differences are larger than the egs_brachy statistical uncertainty.

Also, please repeat the analysis for the central case, too.

[vahx129]

Forgot to add the results for the central case which are here:

Central case (egs_brachy dose - TPS dose) histogram

Central Case (egs_brachy dose - TPS dose)/(relative uncertainty * egs_brachy dose) [denominator is done per voxel]

[mchamberland]

As I'm still thinking about this, here's another suggestion:

How about you simplify as much as possible? Try 1 seed in the center of a phantom in your TPS and in egs_brachy and match the egs++ phantom dimensions exactly to your TPS dose grid. These two dose distributions should match very well.

Then, repeat with the edge case by using a simple water phantom (no egsphant) in egs_brachy (matching your TPS dose grid again) inscribed in a large water box (for scatter purposes).

[vahx129]

Hi again and sorry for the delayed response as I was out of town this last week.

After following what was previously asked for, the results look great for both single seed cases (single central seed and single peripheral seed) even running at 1%/1mm local gamma.

Peripheral Results:

Central Results:

So this makes me wonder if there's something off in my source transformations, superposition settings (which I believe I've fixed based on the previous feedback), and or gamma analysis parameters or expectations. Let me dig into these results a little bit more to start slowly building back up to the full simulation.

[vahx129]

So upon further inspection of my single seed dose distributions, I think saying they looked perfectly fine after just doing a Gamma test was (maybe?) premature. Upon visual inspection of the resulting isodose lines, I see the following:

Coronal Plane MC Dose:

Coronal Plane TPS Dose:

Both Overlaid:

To me, the Monte Carlo distribution looks like it's the one that's more off from what I would expect for a TG-43 superposition simulation. I'm trying to rerun the single seed case with more histories but using1e9 in my original simulation with a single seed was showing relative errors <1% which I don't think fully explain some of those IDLs.

Going to verify the seed model and source distribution values against the library files and then work on dose normalization to see if there's something else I'm messing up here.

[mchamberland]

I agree something looks off on the Monte Carlo dose distribution. It should be symmetric.

Also, can you check if your source is aligned along the same axis as the TPS? It sorta looks like the Monte Carlo model is aligned in the sup-inf direction, but the TPS is ant-post. But maybe something else is going on, too.

[mchamberland]

In your input files, there is a comment that source orientation is not considered in the TPS. Does that mean you're using the anisotropy factor instead of the anisotropy function? It would be more fair to at least align the seeds in the same orientation in both egs_brachy and TPS and to consider its orientation in the TPS.

Also, you should probably disable particle recycling when using superposition mode. In superposition mode, it will only recycle the particles around the currently active source, so even if you rotate the particles, I'm doubtful it will yield much efficiency gains.

And in normal mode, I don't think you gain much efficiency from recycling more than once and from rotating. At least, from my benchmarking, the gains were modest and you're better off setting recycling to 1 and turning off rotation.

[vahx129]

Hi Marc,

I can confirm that the sources are aligned in the same direction in both cases. Please forgive the vague comment in my input (and thanks for reading over it so closely!). What I meant to say is that the actual implanted seeds will have rotations that we do not explicitly model when reconstructing the patient's dose. We manually reconstruct seeds all oriented the same way (long axis aligned in sup/inf directions) centered as closely as possible to the seeds able to visualized in the CT scan.

Example sagittal view of reconstructed seeds in TPS:

Following my TPS coordinate system: +X points towards Patient Left +Y points Posteriorly +Z points Superiorly

TPS TG43 calculation for single seed arbitrarily placed near center of patient CT:

egs_view rendering showing seed long axis is aligned along Z-axis similar to TPS

Also, thank you for the advice regarding particle recycling. I believe for the file I've uploaded into the OneDrive directory called HN6_AWat_5mmXYZ_2cm_Crop_Artifact_1e10_CentralSingleSeed.egsinp that I refer to in this comment and for other files since starting troubleshooting under your guidance have had variance reduction disabled. I verified that all the geometries that I copied into the input were from the CS1 library in egs_brachy and the values matched as well. The CS1 seed also appears as I would expect in egs_view:

I think, in the absence of any other ideas, that I'm going to try to use an EGS_XYZGeometry instead of my egsphant file overridden to water and see if the single seed dose distributions look similar since both the peripheral and central single seed test cases I ran appear to have that weirdly onion shape as opposed to the spherical one we would expect.

[vahx129]

So even in a standard water phantom as opposed to my CT converted to a .egsphant, I'm seeing the same shaped IDL distributions. Do we expect those small points right at the center because of the shape of the seed and maybe less attenuation through the end caps? The overall dose distribution is still more or less spherical with just those peaks near the central axis. The discrepancy in the 5 Gy lines between the TG-43 TPS calculation and the single seed still baffles me a bit. Here's what the 3ddose file looks like in SlicerRT for the single seed in a water phantom XYZ geometry instead of my egsphant files:

The agreement with the TPS single seed calculation is still good overall. I am getting into territory now where selections like the max dose threshold used in the gamma analysis make a big difference in the overall results. Maybe I need to reorient myself and try to figure out what reasonable expectations are for agreement between EGS and TG-43 but I still feel like I'm overlooking something here or possibly getting tripped up somewhere when the clinical implants are being simulated.

I've uploaded the input called WP_Test_1e10.egsinp and its associated 3ddose file to the OneDrive share.

[mchamberland]

I think that shape looks right. Yes, I think the small points are caused by the difference in anisotropy along the source long axis.

How far from the source is the 5 Gy isodose?

[vahx129]

About 3.3 cm

[mchamberland]

The isodoses from your TPS look a little weird, too (for example, the straight edges of the yellow isodose). Is it just an artifact of the isodose display or interpolation?

Along those lines, what's the relationship between the TPS dose grid and the CT grid? Maybe there is a subtlety you're missing here.

[vahx129]

The TPS dose grid uses 1 mm voxels and the egs_brachy simulation uses 0.5 mm voxels. I need to do a little more digging about the way the lines are being represented and will try a software other than SlicerRT to visualize those IDLs to see if that has any underlying implication. I'm also trying to dig into the MIM settings a bit more to try to better understand the CS1 model they're using.

Update: I don't see a ton of things to toggle in MIM but visualizing the same IDLs in the MIM TPS looks more like what I would expect so some of that IDL visualization might just be a rendering issue specific to Slicer. I'm also going to try to use a MATLAB code to run the Gamma analysis so as to also try a different piece of software.

MIM screenshot of single seed in coronal view:

Simultaneously, I'm going to try to create an XYZ geometry that should exactly overlay with the TPS dose grid in terms of dimensions and position and try direct differences to remove the influence of the DTA parameter. Here's another clinical case where the disagreement is high with 5 Gy IDLs shown for both TPS (purple) and egs_brachy (pink). Not only is the general size of the two noticeably different but the Monte Carlo isodose surface also appears to be slightly laterally offset from the TPS one in this axial view. The filled in brown structure is a union of 1 cm spheres at the location of each seed localized in this post implant scan.

But I feel like you're right - there must be some basic subtlety I'm overlooking here but I'm surprised it only appears to happen for some clinical cases and not others. The case from the screenshot in this post utilized 28 seeds for what it's worth which is definitely on the higher side for these cases. Though I'm not sure if I can generally state that more seeds leads to more discrepancy without a bit more digging. .

[mchamberland]

Any update?

I’m in the process of submitting a proposal to my IRB for a retrospective Monte Carlo study of our prostate seed implant cases. I’m curious to see if I encounter similar issues in that context and with VariSeed as the TPS.

[vahx129]

Hi again Marc,

Here's a (very delayed, sorry!) update on what I've been up to so far. I've been resimulating the TG-43 comparisons using the suggestions you provided previously and trying to document the results. Based on what I'm seeing so far, I'm honestly not sure if seed position relative to center/periphery is actually what's causing the discrepancies. Some comparisons look fantastic and others look really bad and I'm honestly not totally sure what stands out about the ones that look bad. Here's a summary table of where I'm at so far.

To quickly recap, the TG-43 simulations are run in egs_brachy using voxels as close to 0.5 mm x 0.5 mm x 0.5 mm that I can get away with. If the CT voxels are larger than 0.5 mm in any specific direction, I don't generally resample and just accept whatever that value is. In most cases, the CT resolution will be 0.5 to 0.6 mm in all directions. I use the CT images to generate a .egsphant file for each geometry and use a ramp that makes all voxels into water. I will usually crop the CT image using ctcreate to only include the slices which have the seeds and an extra couple of cm superior and inferior to save on simulation time. I then convert each RT-Plan file into egs_brachy compatible inputs and use a normalization based on the CS-1 source from Carlton and validated by my own independent TG-43 simulation. After the simulations are completed, I import the 3ddose file into SlicerRT alongside the RT-Dose file from our TPS (MIM) and run the comparisons. I am using Global Gamma tests that compare values as a percentage of 60 Gy (which is our standard Rx dose for these cases) and keep a 10% dose threshold to ignore all voxels lower than 6 Gy.

If there are any specific simulation files or anonymized doses I can provide to help, please let me know. I can also try to look at the uncertainties on the simulations themselves (I'm now running 10x more histories than I was previously).

[vahx129]

I believe that I may have identified the error in my methodology but I can't say for sure just yet. I'm going to run a handful of the simulations with the worst results using the DICOM CT and Dose data resampled to 0.5 mm cubic voxels. At least one simulation using properly matched voxel sizes appears to be agreeing much better already. Will reach out again when I have some more of the results ready.

[vahx129]

This appears to have been a total user error on my part. My ctcreate all water ramp accidentally retained a line from another file, which allowed densities other than 1 g/cc to be assigned to CT voxels. So, in cases where the seeds were almost entirely surrounded by tissue, results were good, whereas seeds near the edge of the patient's anatomy and/or bones looked bad because of the change in density despite the material being assigned to water.

Sorry this didn't occur to me sooner. I hadn't gotten dosxyz_show to work on my cluster but when I took a look on my personal PC, this became more apparent. Since my cluster is down at the moment, I can't say with 100% certainty that this is not a bug but I'm about 99.8% sure this is not actually a bug.

[mchamberland]

No problem at all. It's always those little things that are hard to track down.