| Geant4 Cross Reference |
1 #!/usr/bin/env python
2 # coding: utf-8
3
4 # To run, change the links and definitions in this script and
5 # type in the terminal containing the file
6 # $python3 molecularDNAsurvival.py
7
8 # Authors:
9 # K. Chatzipapas (*), D. Sakata, H. N. Tran, M. Dordevic, S. Incerti et al.
10 # (*) contact: k.chatzipapas@yahoo.com
11
12 # Authors' note : This is a very early version of this calculation technique.
13 # Needs to optimized for each different application. An updated
14 # version will be released in the future.
15
16
17 # Define filenames
18 outputFile = "molecularDNAsurvival.txt"
19 iRootFile = "../molecular-dna.root"
20 print("\nInput file:",iRootFile,"\n")
21
22 # Define cell parameters
23 r3 = 7100*1e-09 * 2500*1e-09 * 7100*1e-09 # a * b * c / Chromosome size, as defined in the mac file, but in meters. If sphere, a=b=c
24 NBP = 6405886128 # bp / Length of the DNA chain in Mbp
25 mass = 997 * 4 * 3.141592 * r3 / 3 # waterDensity * 4/3 * pi * r3 in kg
26
27 # Integration parameters
28 t0 = 0.0 # Starting time point (dimensionless)
29 t1 = 336 # Final time point (dimensionless)
30 dt = 0.001 # Intergration time step (dimensionless)
31 numpoints = int( (t1 - t0) / dt )
32
33 # Survival model Parameters
34 # THESE PARAMETERS HAVE TO BE DEFINED FOR PROPER RESULTS/ EXPERIMENTAL PARAMETERS,
35 # AS DESCRIBED IN (Stewart RD. Radiat Res. 2001 https://pubmed.ncbi.nlm.nih.gov/11554848/ )
36 maxDose = 9
37 DR1 = 720 #Gy/h SF-dose
38 bp = 1 #Gbp
39 #DR2 = 600. #Gy/h SF-time
40 #DR3 = 12000. #Gy/h FAR
41
42 # usual value of gamma 0.25 (changes the end value of the curve)
43 gamma = 0.19 #other published values: 0.189328
44 # Lamb1 the end of the curve Lamb2 is connected with Eta
45 Lamb1 = 0.10 #other published values: 0.0125959 # 0.671 h-1 or 2.77 h-1 (15 min repair half-time)
46 Lamb2 = 1 #other published values: 33062.9 # 0.0439 h-1 (15.8-h repair half-time) or 0.0616 h-1 (11.3 h repair half-time)
47 #
48 Beta1 = 0 #other published values: 0.0193207 # 0.00152
49 Beta2 = 0
50 # Defines the curvature as well as the end value of the curve
51 Eta = 0.00005 #other published values: #7.50595e-06 # 1.18e-04 h-1 or 0.00042 h-1
52 # name of the cell (used in the output)
53 cell = "test"
54
55 # ONLY ADVANCED USERS AFTER THIS POINT
56 #####################################################################################################################
57 import ROOT as root
58 import numpy as np
59 import matplotlib.pyplot as plt
60 from scipy.integrate import ode
61 import threading
62
63 # Definition of the survival model
64 def SF_system( T_, DR_, Y_, Sig1_, Sig2_, gamma_, Lamb1_, Lamb2_, Beta1_, Beta2_, Eta_ ):
65
66 global T
67 global DR
68 global Y
69 global Sig1
70 global Sig2
71 global gamma
72 global Lamb1
73 global Lamb2
74 global Beta1
75 global Beta2
76 global Eta
77
78 T = T_
79 DR = DR_
80 Y = Y_
81 Sig1 = Sig1_
82 Sig2 = Sig2_
83 gamma = gamma_
84 Lamb1 = Lamb1_
85 Lamb2 = Lamb2_
86 Beta1 = Beta1_
87 Beta2 = Beta2_
88 Eta = Eta_
89
90
91 return T, DR, Y, Sig1, Sig2, gamma, Lamb1, Lamb2, Beta1, Beta2, Eta
92
93 def SF_function( t, x, p ):
94 if(t<=T):
95 #print("t<=T")
96 dx[0] = DR*Y*Sig1 - Lamb1*x[0] - Eta*x[0]*(x[0]+x[1])
97 dx[1] = DR*Y*Sig2 - Lamb2*x[1] - Eta*x[1]*(x[0]+x[1])
98 dx[2] = Beta1*Lamb1*x[0] + Beta2*Lamb2*x[1]+gamma*Eta*(x[0]+x[1])*(x[0]+x[1])
99 dx[3] = (1.-Beta1)*Lamb1*x[0] + (1.-Beta2)*Lamb2*x[1]+(1.-gamma)*Eta*(x[0]+x[1])*(x[0]+x[1])
100 else:
101 dx[0] = - Lamb1*x[0] - Eta*x[0]*(x[0]+x[1])
102 dx[1] = - Lamb2*x[1] - Eta*x[1]*(x[0]+x[1])
103 dx[2] = Beta1*Lamb1*x[0] + Beta2*Lamb2*x[1]+gamma*Eta*(x[0]+x[1])*(x[0]+x[1])
104 dx[3] = (1.-Beta1)*Lamb1*x[0] + (1.-Beta2)*Lamb2*x[1]+(1.-gamma)*Eta*(x[0]+x[1])*(x[0]+x[1])
105 return dx
106
107
108 DSBBPID = []
109 repDSB = 0
110 irrDSB = 0
111 acc_edep = 0
112 dose = 0
113 i = 0
114 eVtoJ = 1.60218e-19
115
116
117 # Analyze root files
118 f = root.TFile(iRootFile)
119
120 fTree = f.Get("tuples/damage")
121 for e in fTree:
122 if (e.TypeClassification=="DSB" or e.TypeClassification=="DSB+" or e.TypeClassification=="DSB++"):
123 DSBBPID.append((e.Event, e.BasePair))
124 i += 1
125
126 gTree = f.Get("tuples/classification")
127 for e in gTree:
128 repDSB += e.DSB
129 irrDSB += e.DSBp + 2*e.DSBpp
130
131 ffTree = f.Get("tuples/chromosome_hits")
132 for ev in ffTree:
133 acc_edep += (ev.e_chromosome_kev + ev.e_dna_kev) *1e3 # eV
134
135
136 # Calculate the total number of DSBs
137 totalDSB = repDSB + irrDSB
138 # Calculate the absorbed dose
139 dose = acc_edep * eVtoJ / mass # Dose in Gy
140
141 #print ("r3=", r3)
142 #print("Dose :", dose, "\nTotal DSB/Gy :", totalDSB/dose, "\nRepairable DSBs/Gy :", repDSB/dose, "\nIrrepairable DSBs/Gy :", irrDSB/dose)
143
144 NBP_of_cell = 4682000000 # This is a constant of the model, as defined by Stewart in the TLK (Stewart RD. Radiat Res. 2001 https://pubmed.ncbi.nlm.nih.gov/11554848/ )
145 NBPcell = NBP/NBP_of_cell
146 totalDSBperGyCell = totalDSB/NBPcell/dose
147 #totalDSBperGyCell = totalDSB/dose
148 print ("Dose :", dose, "\nTotal DSB/Gy/Cell :", totalDSBperGyCell)
149
150 kFactor = 1 # This factor is used to normalize irrepairable damage to lower values than 1.
151 S1 = repDSB/NBPcell/dose/kFactor
152 S2 = irrDSB/NBPcell/dose/kFactor
153 print("Repairable DSBs/Gy/cell :", S1, "\nIrrepairable DSBs/Gy/cell :", S2, "\n\n")
154
155
156 # In[ ]:
157
158
159 T = 0
160 DR = 0
161 Y = 0
162 kozi = False
163 sum = [0]*(maxDose+1)
164 outputfile = [""] * (maxDose+1)
165
166 # Create the time samples for the output of the ODE solver.
167 # I use a large number of points, only because I want to make
168 # a plot of the solution that looks nice.
169 #t = [(t1-t0) * float(i) / (numpoints - 1) for i in range(numpoints)]
170 t = np.linspace(t0, t1, numpoints)
171
172 dx = np.zeros(4)
173 x = np.zeros(4)
174 dsbs = 0.0
175 Sig1 = S1
176 Sig2 = S2
177
178 filename = "txt/molecularDNAsurvival_"+cell+".dat"
179 outputfile1 = open(filename,"w")
180 outputfile1.write("Dose \t Survival Fraction\n")
181
182 for j in range(maxDose+1):
183 #j=1
184 sol = []
185 time = []
186 dose = j
187 expEndTime = dose/DR1 #h
188
189 #////double StopTime = (expEndTime+ n/min(Lamb1,Lamb2));
190
191 arguments = SF_system( expEndTime, DR1,bp,Sig1,Sig2,gamma,Lamb1,Lamb2,Beta1,Beta2,Eta )
192 p = [ T, DR, Y, Sig1, Sig2, gamma, Lamb1, Lamb2, Beta1, Beta2, Eta ]
193
194 Stepper = "dopri5" # dopri5,vode,dop853
195
196 if (j==0):
197 solver = ode(SF_function).set_integrator(Stepper, nsteps=1)
198 print("#*#*#*#*#*# Ignore this warning #*#*#*#*#*#")
199 else:
200 solver = ode(SF_function).set_integrator(Stepper, nsteps=numpoints)
201 solver.set_initial_value( x, t0 ).set_f_params(p)
202 #print(solver.t,solver.y)
203
204 # Repeatedly call the `integrate` method to advance the
205 # solution to time t+dt, and save the solution in solver.
206 while solver.t < t1: # and solver.successful():
207 solver.integrate(solver.t+dt)
208
209 sum[j] = solver.y[2]
210 print("Survival probability of cell ("+cell+") after",j,"Gy is", np.exp(-sum[j]) )
211 outputfile1.write(str(dose) +"\t"+ str(np.exp(-(solver.y[2]))) +"\n")
212
213 c = np.linspace(0, maxDose, num=maxDose+1)
214 sum=np.array(sum)
215 outputfile1.close()
216
217 print("\nCalculation concluded. Check the final image and close it to end the program.!")
218
219 plt.yscale('log')
220 plt.xlim([0, 8])
221 plt.ylim([1e-2, 1])
222 plt.plot(c,np.exp(-sum[:]))
223 plt.show()
224
225 print("\nThank you for using molecularDNAsurvival.!\n")
226