[Run issue]: Step-by-Step in TsNucleus possible Scorer/ Physics and chemistry problem

[ShahinGHD]

Contact Details

[email protected]

What was the issue?

Hi everyone,
This issue needs to be read carefully. I apologize for the long text.
it's been a while since I have tried to reproduce the results from the following TOPAS-nBio paper:

Impact of DNA Geometry and Scoring

I am using a single threshold method for direct damages, which is 11.75 ev for the Sphere model and 17.5 ev for the half-cylinder model,
I used the probability of scavenging in backbones of 0.0585 for the half cylinder, and 0.25 for the sphere model as well for indirects.
I believe I have the exact beam setup with 1MeV protons (LET ~ 29 KeV/um) shooting from the edge of the nucleus with 10000 histories.
Please note that instead of including the "DNAScorerStepByStep.txt" file, I entered the values in the attached input file to show how I have done it.
I did several runs and rotated the nucleus for each run.

I encountered several issues with these models in physics and chemistry.

Problem 1 : Direct Damage recording

I tried to turn off the chemistry and only use "TsEmDNAPhysics" to reproduce the "direct" BD/Gy/Gbp and SB/Gy/Gbp from Figure 3 in the paper.
Figure 3 shows:

half-cylinder model "direct" SB/Gy/Gbp = 50-60 ||==> my simulation = 26 (column 5 in the output file)
half-cylinder model "direct" BD/Gy/Gbp = 10-13 ||==> my simulation = 15.9 (column 9 in the output file)

The output Excel file is attached below: (half-cylinder model, Only Direct damage recorded - Chemistry has been turned off)
1 MeV Proton beam uniform field at 4.65 um from nucleus center as mentioned in the paper.

half C direct only 10000H.xlsx

Problem 2 : Simulation with chemistry transport

Since the full direct and indirect simulations are slow, I had 100 histories per run when chemistry was on but I did multiple runs that took 20h each.
The Chemical High end was set to 1 ns. sampling at 0.5 PS as mentioned in the paper.

For the Sphere model it was said the ratio of BD/SB for "indirect" should be ~2.5-3.

using the "TsEmChemistry" module and "TOPASDefaultReactions.txt" file, I could never hit this ratio, it was always less which was (~ 1.8).
I tried to use the "Extended Chemistry" and "TOPASChemistry_Extended.txt". I could obtain the ratio of "2.4" with 100 histories in 2 or 3 runs, but I could never reproduce Figure 7 for both the sphere and half-cylinder models.

for example in Figure 7, the TOTAL SSB/Gy/Gbp is: (column 4 in the output) (direct+indirect)
half-cylinder ===> paper =~ 95,||==> my simulation ~ 31
Sphere model ==> paper = ~70,||==> my simulation ~ 18

it seems the indirect damages recorded are underestimated here compared to the paper by a large amount.

Problem 3 : Comparison between SBS "Direct" DNA damages and DBSCAN using linear probability

Comparing Step by step direct DNA damages and the DBSCAN, both with "linear probability" which is 5ev to 37.5 ev shows huge differences. 1Mev Protons, 1000 History. same beam setup. Here is the averaged number of SSB and DSB recorded among all events in both methods compared:

DNA MODEL : Sphere

Step-by-step (SSB) = 14.9 ==> column 17 in the output file.
Step-by-step (DSB) = 0.83 ==> column 10 in the output file

DBSCAN (SSB) mean = 323
DBSCAN (DSB) mean = 35
DBSCAN (complex DSB) = 28

DNA MODEL: Half-Cylinder

Step-by-step (SSB) = ==> 37.54 column 17 in the output file.
Step-by-step (DSB) = ==> 5.73 column 10 in the output file

DBSCAN (SSB) mean = 330
DBSCAN (DSB) mean = 36
DBSCAN (complex DSB) = 29

Can you please be kind and look at the input file and also confirm that the chemistry and physics modules as well as the DNA models are working correctly for the TSNucleus?
Any information would be appreciated.

PS: Visualization is also not available and QT crashes. also, it does not even show outlines when options are set in the visualization sections.

Please and thank you,
Shahin

TOPAS-nBio category

Chemistry

Version

3.0

What Operative System (OS) are you running on?

Windows (WSL)

Relevant log output

#############
# GEOMETRY #
#############


# World
d:Ge/World/HLX                          = 14 um
d:Ge/World/HLY                          = 14 um
d:Ge/World/HLZ                          = 14 um
b:Ge/World/Invisible                    = "true"
s:Ge/World/Material                     = "G4_WATER"

# Nucleus
s:Ge/Nucleus/Type                       = "TsNucleus"
s:Ge/Nucleus/Parent                     = "World"
s:Ge/Nucleus/Material                   = "G4_WATER"
s:Ge/Nucleus/Color                      = "white"
s:Ge/Nucleus/DNAModel					= "HalfCylinder"
d:Ge/Nucleus/NucleusRadius              = 4.65 um
i:Ge/Nucleus/HilbertCurveLayer          = 4
i:Ge/Nucleus/HilbertCurve3DRepeat       = 30
s:Ge/Nucleus/HilbertCurveFileName       = "supportFiles/HilbertPoints_iteration1.dat"
s:Ge/Nucleus/Voxel/Material             = "G4_WATER"
b:Ge/Nucleus/ShowDNAVoxels              = "F"
s:Ge/Nucleus/ChromatinFiber/Material    = "G4_WATER"
b:Ge/Nucleus/ShowChromatinCylinders     = "F"
s:Ge/Nucleus/Histone/Material           = "G4_HistoneMaterial"
b:Ge/Nucleus/OnlyBuildOneHistone        = "false"
b:Ge/Nucleus/HistoneAsScavenger         = "true"
b:Ge/Nucleus/AddBases                   = "true"
s:Ge/Nucleus/Base1/Material             = "G4_BaseMaterial"
s:Ge/Nucleus/Base2/Material             = "G4_BaseMaterial"
b:Ge/Nucleus/AddBackbones               = "true"
s:Ge/Nucleus/Backbone1/Material         = "G4_BackboneMaterial"
s:Ge/Nucleus/Backbone2/Material         = "G4_BackboneMaterial"
b:Ge/Nucleus/AddHydrationShell          = "true"
s:Ge/Nucleus/HydrationShell1/Material   = "G4_WATER_MODIFIED"
s:Ge/Nucleus/HydrationShell2/Material   = "G4_WATER_MODIFIED"
# d:Ge/Nucleus/fHydrationShellThickness = 0.16 nm
#
# # -- Materials
s:Ma/G4_HistoneMaterial/CloneFromMaterial   = "G4_WATER"
d:Ma/G4_HistoneMaterial/CloneWithDensity    = 1.407 g/cm3
s:Ma/G4_BaseMaterial/CloneFromMaterial      = "G4_WATER"
d:Ma/G4_BaseMaterial/CloneWithDensity       = 1.0 g/cm3
s:Ma/G4_BackboneMaterial/CloneFromMaterial  = "G4_WATER"
d:Ma/G4_BackboneMaterial/CloneWithDensity   = 1.0 g/cm3
s:Ma/G4_WATER_MODIFIED/CloneFromMaterial    = "G4_WATER"
d:Ma/G4_WATER_MODIFIED/CloneWithDensity     = 1.0 g/cm3


b:Ge/Nucleus/RotateNucleusForEachRun				= "True"




# SOURCE

s:So/ParticleGun/Type                               = "Beam"
s:So/ParticleGun/Component                          = "BeamPosition"

s:So/ParticleGun/BeamParticle                       = "Proton"				 
s:So/ParticleGun/BeamPositionDistribution			= "Flat"					
s:So/ParticleGun/BeamPositionCutoffShape			= "Ellipse"
d:So/ParticleGun/BeamPositionCutoffX 				= Ge/Nucleus/NucleusRadius um
d:So/ParticleGun/BeamPositionCutoffY 				= Ge/Nucleus/NucleusRadius um
s:So/ParticleGun/BeamAngularDistribution            = "none"
d:So/ParticleGun/BeamEnergy                         = 1 MeV
i:So/ParticleGun/NumberOfHistoriesInRun             = 10000

d:Ge/BeamPosition/TransZ = 4.65 um


# PHYSICS

sv:Ph/Default/Modules       =  1 "TsEmDNAPhysics" #"TsEmDNAChemistry" # turned off chemistry
s:Ph/Default/Electron/SetElasticScatteringModel = "CPA100"

# This section was commented out for direct DNA break measurements
# includeFile						=supportFiles/TOPASDefaultReactions.txt
# s:Ch/ChemistryName = "TOPASChemistry"
# b:Ch/TOPASChemistry/ChemicalStageTransportActive = "F"
# d:Ch/TOPASChemistry/ChemicalStageTimeEnd = 1 ns
# dv:Ch/TOPASChemistry/ChemicalStageTimeStepsHighEdges = 1 1 ns
# dv:Ch/TOPASChemistry/ChemicalStageTimeStepsResolutions = 1 0.5 ps



# SCORER

# includeFile					= supportFiles/DNAScorerStepByStep.txt  # commented out
s:Sc/DNADamageScorer/Quantity				= "DNADamageNucleusStepByStep"
s:Sc/DNADamageScorer/Component 				= "Nucleus"
i:Sc/DNADamageScorer/NumberOfHistoriesInRun = So/ParticleGun/NumberOfHistoriesInRun
d:Sc/DNADamageScorer/ScoringRadius			= Ge/Nucleus/NucleusRadius um
i:Sc/DNADamageScorer/BasePairPositionAtGeometricHierarchy = 0
sv:Sc/DNADamageScorer/Strand1MaterialNames  = 4 "G4_WATER_MODIFIED" "G4_BackboneMaterial" "G4_BaseMaterial" "G4_HistoneMaterial"
sv:Sc/DNADamageScorer/Strand2MaterialNames  = 3 "G4_WATER_MODIFIED" "G4_BackboneMaterial" "G4_BaseMaterial"
s:Sc/DNADamageScorer/PrimaryParticle		= So/ParticleGun/BeamParticle
d:Sc/DNADamageScorer/MeanEnergy				= So/ParticleGun/BeamEnergy MeV
# d:Sc/DNADamageScorer/StopTrackingAtDose		= 0.5 gray
d:Sc/DNADamageScorer/ImageXMinPosition		= 0 um - Ge/Nucleus/NucleusRadius
d:Sc/DNADamageScorer/ImageXMaxPosition		= Ge/Nucleus/NucleusRadius um
d:Sc/DNADamageScorer/ImageYMinPosition		= 0 um - Ge/Nucleus/NucleusRadius
d:Sc/DNADamageScorer/ImageYMaxPosition		= Ge/Nucleus/NucleusRadius um
d:Sc/DNADamageScorer/ImageZMinPosition		= 0 um - Ge/Nucleus/NucleusRadius
d:Sc/DNADamageScorer/ImageZMaxPosition		= Ge/Nucleus/NucleusRadius um


# Setting up the output
b:Sc/DNADamageScorer/PropagateToChildren						= "true"
s:Sc/DNADamageScorer/OutputFile									= "DNADamage"
s:Sc/DNADamageScorer/OutputType									= "ascii"
s:Sc/DNADamageScorer/IfOutputFileAlreadyExists					= "overwrite"

# Options for direct damage (recommended values for spherical DNA model)
d:Sc/DNADamageScorer/DirectDamageThreshold						= 17.5 eV    # For half-cylindrical DNA model, the recommended value is 17.5 eV
b:Sc/DNADamageScorer/UseLinearProbabilityForDirectDamage		= "false"
#d:Sc/DNADamageScorer/LowerLimitForLinearProbabilityFunction	= 5 eV
#d:Sc/DNADamageScorer/UpperLimitForLinearProbabilityFunction	= 37.5 eV

# Options for quasi-direct damage
u:Sc/DNADamageScorer/ProbabilityOfChargeTransferFromHydrationShellToBackbone = 0.33333

# Options for indirect damage (recommended values for spherical DNA model)
# b:Sc/DNADamageScorer/AlwaysScavengeSpeciesInDNAComponents		= "F"
# # u:Sc/DNADamageScorer/ProbabilityForOHToCauseDamage = 0.4
# u:Sc/DNADamageScorer/ProbabilityOfScavengingInBackbone			= 0.0585		  # For half-cylindrical DNA model, the recommended value is 0.0585
# u:Sc/DNADamageScorer/ProbabilityOfScavengingInBase				= 1.0
# u:Sc/DNADamageScorer/ProbabilityOfIndirectDamageToBackbone		= 0.55		  # If always scavenging in DNA compoennts, the recommended value is 0.4
# u:Sc/DNADamageScorer/ProbabilityOfIndirectDamageToBase			= 1.0		  # If always scavenging in DNA compoennts, the recommended value is 0.4
# b:Sc/DNADamageScorer/ScavengeInHistones							= "true"

# Options for damage classification
i:Sc/DNADamageScorer/MaximumBasePairDistanceToConsiderDSB		= 10
b:Sc/DNADamageScorer/ExcludeShortFragments						= "false"
# i:Sc/DNADamageScorer/LowerThresholdForFragmentDetection		= 0
# i:Sc/DNADamageScorer/UpperThresholdForFragmentDetection		= 3E8

# Options for foci scoring
b:Sc/DNADamageScorer/ScoreNumberOfFoci							= "true"
dv:Sc/DNADamageScorer/FociSizes									= 2 150 500 nm # Maximum of 5 sizes
b:Sc/DNADamageScorer/Get3DFociImage								= "true"
b:Sc/DNADamageScorer/Get2DFociImages							= "true"
sv:Sc/DNADamageScorer/2DFociImagePlanes							= 3 "X" "Y" "Z"
s:Sc/DNADamageScorer/MicroscopePSFShape							= "Gaussian"
d:Sc/DNADamageScorer/MicroscopePSFWidth							= 400 nm
dv:Sc/DNADamageScorer/Foci2DImageResolutions					= 2 25 100 nm
d:Sc/DNADamageScorer/Foci3DImageResolution						= 400 nm

# Options for output extension
b:Sc/DNADamageScorer/WriteCSVOutputWithAllDamageSpecification	= "true"
b:Sc/DNADamageScorer/ScoreDirectDamage							= "true"
b:Sc/DNADamageScorer/ScoreIndirectDamage						= "true"
b:Sc/DNADamageScorer/ScoreQuasiDirectDamage						= "true"
b:Sc/DNADamageScorer/ScoreBaseDamages							= "true"
b:Sc/DNADamageScorer/ScoreBackboneDamages						= "true"
b:Sc/DNADamageScorer/BreakDownOutputPerDamageOrigin				= "true"

# Options for SDD output specification
d:Sc/DNADamageScorer/DosePerExposure							= 1 Gy
b:Sc/DNADamageScorer/IncludeDSBDamageSitesOnlyInSDD				= "false"
b:Sc/DNADamageScorer/WriteMinimalSDDOutput						= "false"

# Data for SDD header
s:Sc/DNADamageScorer/AuthorName							   		= "Name"
s:Sc/DNADamageScorer/SimulationDetails					   		= "Nucleus simulation"
s:Sc/DNADamageScorer/SourceDetails					   	   		= "Monoenergetic source"
i:Sc/DNADamageScorer/SourceType					   		   		= 1
s:Sc/DNADamageScorer/EnergyDistribution					   		= "M, 0"
s:Sc/DNADamageScorer/IrradiationTarget					   		= "Nucleus"
s:Sc/DNADamageScorer/CellCycleStage					       		= "0"
s:Sc/DNADamageScorer/DNAStructure					           	= "0, 1"
i:Sc/DNADamageScorer/InVitroOrInVivo						   	= 0
s:Sc/DNADamageScorer/ProliferationStatus					   	= "1"
s:Sc/DNADamageScorer/Microenvironment					       	= "20, 0.01"
d:Sc/DNADamageScorer/Time					  				   	= 0 s
s:Sc/DNADamageScorer/AdditionalInfo					       		= ""





# Geometry file for Nucleus scorer
s:Sc/DNADamageScorer/GeometryInfo			= "supportFiles/Geoinfo.txt"
s:Sc/DNADamageScorer/FileVoxelCopyNumbersInsideNucleus = "supportFiles/CopyNoTable.txt"
s:Sc/DNADamageScorer/MapChromosomeVoxelsFile = "supportFiles/signedCHVoxel.txt"

i:Ts/ShowHistoryaCountAtInterval                       = 1



# VISUALIZATION

# b:Gr/Enable = "T"
# s:Gr/View/Type			= "OpenGL"
# s:Gr/View/ColorBy       = "ParticleType"
# # s:Gr/View/ColorBy       = "charge"
# # sv:Gr/View/ColorByChargeColors = 3 "yellow" "aqua" "blue"
# # # sv:Gr/View/ColorByParticleTypeNames = 5 "e-" "gamma" "proton" "alpha" "neutron"
# # # sv:Gr/View/ColorByParticleTypeColors = 5 "green" "yellow" "blue" "red" "brown"
# sv:Gr/view/ColorByParticleTypeNames = 12 "alpha" "e-" "gamma" "Hydroxil" "SolvatedElectron" "HydrogenPeroxide" "Hydrogen" "DiHydrogen" "Oxygen" "SuperoxideAnion" "Hydroperoxy" "Dioxidanide"
# sv:Gr/view/ColorByParticleTypeColors = 12 "yellow" "red" "green" "violet" "magenta" "grass" "blue" "aqua" "white" "purple" "lightblue" "blue"
# b:Ts/UseQt              = "True"
# i:Gr/ShowOnlyOutlineIfVoxelCountExceeds = 10
# # i:Gr/SwitchOGLtoOGLIifVoxelCountExceeds = 1


i:Ts/Seed = 4124
i:Ts/NumberOfThreads = 0

Code of Conduct

• I agree to follow this project's Code of Conduct

Attestation

• I testify that I have fully read the bug reporting guideline beforehand.

[ShahinGHD]

hello everyone,
can someone please confirm if the step by step method is not working properly?
as I discussed above, it is not possible to reproduce the mentioned paper results,

Thanks
Shahin