| Geant4 Cross Reference |
1 #!/usr/bin/env python
2 # coding: utf-8
3
4 # To run, change the links and definitions in
5 # type in the terminal containing the file
6 # $python3 molecularDNArepair.py
7
8 # Authors:
9 # K. Chatzipapas (*), D. Sakata, H. N. Tran, M
10 # (*) contact: k.chatzipapas@yahoo.com
11
12 # Define filenames
13 outputFile = "molecularDNArepair.txt"
14 iRootFile = "../molecular-dna.root"
15 print("\nInput file:",iRootFile,"\n")
16
17 # Define cell parameters
18 r3 = 7100*1e-09 * 2500*1e-09 * 7100*1e-09 # a
19 NBP = 6405886128 # bp / Length of the DNA chai
20 mass = 997 * 4 * 3.141592 * r3 / 3 # water
21
22 # Integration parameters
23 t0 = 0.0; # Starting time point (dimensionle
24 t1 = 25; # Final time point (dimensionless)
25 dt = 0.001 #0.001; # Intergration time step (d
26 numpoints = int( (t1 - t0) / dt )
27
28
29 # ONLY ADVANCED USERS AFTER THIS POINT
30 # Repair model parameters may be find at line
31 ##############################################
32 import ROOT as root
33 import numpy as np
34 import matplotlib.pyplot as plt
35 from scipy.integrate import odeint as od
36
37 # Initialization
38 DSBBPID = []
39 i = 0
40 repDSB = 0
41 irrDSB = 0
42 acc_edep = 0
43 dose = 0
44 eVtoJ = 1.60218e-19
45
46
47 f = root.TFile(iRootFile)
48
49 fTree = f.Get("tuples/damage")
50 for e in fTree:
51 if (e.TypeClassification=="DSB" or e.TypeC
52 DSBBPID.append((e.Event, e.BasePair))
53 i += 1
54
55 gTree = f.Get("tuples/classification")
56 for e in gTree:
57 repDSB += e.DSB
58 irrDSB += e.DSBp + e.DSBpp # 2*eDSBpp c
59
60 ffTree = f.Get("tuples/chromosome_hits")
61 for ev in ffTree:
62 acc_edep += (ev.e_chromosome_kev + ev.e_dn
63
64
65
66 # Calculate the total number of DSBs
67 totalDSB = repDSB + irrDSB
68
69 # Calculate the absorbed dose
70 dose = acc_edep * eVtoJ / mass # Dose in Gy
71
72 #print ("r3=", r3)
73 #print("Dose :", dose, "\ni :", i, "\nTotal DS
74
75 NBP_of_cell = 5000000000 # This is a constan
76 NBPcell = NBP/NBP_of_cell
77 totalDSBperGyCell = totalDSB/NBPcell/dose
78 print ("Dose :", dose, "\nTotal DSB/Gy/Cell :"
79 S1 = repDSB/NBPcell/dose
80 S2 = irrDSB/NBPcell/dose
81 print("Repairable DSBs/Gy/cell :", S1, "\nIrre
82
83
84 # Model parameters
85 fDz = 1 # Dose is set to 1 Gy, because fal
86 falpha = totalDSBperGyCell #27.5 #totalDSB #
87 fNirrep = irrDSB/totalDSB # irrepDSB # racirre
88 if (fNirrep == 0):
89 fNirrep = 0.01 # For not printing error
90 print("Ratio of irrepairable to repairable dam
91
92 # Initial conditions for ODE
93 NbEquat = 29
94 Y = [0] * NbEquat
95 #print (Y[1])
96 Y[0] = falpha * fDz
97 YP = [0] * NbEquat
98 print("\nY[0] used in the model :", Y[0])
99
100 # Create the time samples for the output of th
101 # I use a large number of points, only because
102 # a plot of the solution that looks nice.
103 t = [(t1-t0) * float(i) / (numpoints - 1) for
104
105
106 # DIMENSIONAL REACTION RATES DEFINITON (check
107 #
108 #------------NHEJ--------------
109 K1 = 11.052 #other:10.02 #11.052
110 Kmin1 = 6.6*1e-01 #6.59999*1e-04 # h-1
111 K2 = 5.82*1e+05 #18.8305*(1.08517-np.exp(-21
112 #print (K2)
113 Kmin2 = 5.26*1e-01 #h-1
114 K3 = 1.86 # h-1
115 K4 = 1.38*1e+06 # M-1*h-1
116 Kmin4 = 3.86*1e-04 # h-1
117 K5 = 15.24 # M-1*h-1
118 Kmin5 = 8.28 # h-1
119 K6 = 18.06 # M-1*h-1
120 Kmin6 = 1.33 # h-1
121 K7 = 2.73*1e+05 # M-1*h-1
122 Kmin7 = 3.2 # h-1
123 #K8 is a normalization factor and sometimes it
124 normK8 = 0.35
125 K8 = normK8*5.52*1e-01 # h-1
126 K9 = 1.66*1e-01 # h-1
127 K10 = (1.93*1e-07)/fNirrep # M
128 K11 = 7.50*1e-02 # h-1
129 K12 = 11.1 # h-1
130 #
131 #------------HR--------------
132 P1 = 1.75*1e+03 # M-1*h-1
133 Pmin1 = 1.33*1e-04 # h-1
134 P2 = 0.39192 # h-1
135 Pmin2 = 2.7605512*1e+02 # h-1
136 P3 = 1.37*1e+04 # M-1*h-1
137 Pmin3 = 2.34 # h-1
138 P4 = 5.52*1e-2 #3.588*1e-02 # h-1
139 P5 = 1.20*1e+05 # M-1*h-1
140 Pmin5 = 8.82*1e-05 # h-1
141 P6 = 1.87*1e+05 #1.54368*1e+06 # M-1*h
142 Pmin6 = 1.55*1e-03 # h-1
143 P7 = 21.36 #1.4904 # h-1
144 P8 = 1.20*1e+04 # M-1*h-1
145 Pmin8 = 2.49*1e-04 # h-1
146 P9 = 4.88*1e-1 #1.104 # h-1
147 P10 = 7.20*1e-03 # h-1
148 P11 = 6.06*1e-04 # h-1
149 P12 = 2.76*1e-01 # h-1
150 #
151 #------------SSA--------------
152 Q1 = 7.8*1e+03 #1.9941*1e+05 # M-1*
153 Qmin1 = 1.71*1e-04 # h-1
154 Q2 = 3*1e+04 #4.8052*1e+04 # M-1*
155 Q3 = 6*1e+03 # M-1*h-1
156 Qmin3 = 6.06*1e-04 # h-1
157 Q4 = 1.66*1e-06 #1.62*1e-03 # h-1
158 Q5 = 8.40*1e+04 # M-1*h-1
159 Qmin5 = 4.75*1e-04 # h-1
160 Q6 = 11.58 # h-1
161 #
162 #-------alt-NHEJ (MMEJ)--------
163 R1 = 2.39*1e+03 # M-1*h-1
164 Rmin1 = 12.63 # h-1
165 R2 = 4.07*1e+04 # M-1*h-1
166 R3 = 9.82 # h-1
167 R4 = 1.47*1e+05 # M-1*h-1
168 Rmin4 = 2.72 # h-1
169 R5 = 1.65*1e-01 # h-1
170 #
171 # Scalling rate XX1
172 XX1 = 9.19*1e-07 # M
173 XX3 = 2.3 *1e-12 # M
174 #
175 # DIMENSIONLESS REACTION RATES
176 #
177 #------------NHEJ--------------
178 k1 = K1*XX1/K8
179 kmin1 = Kmin1/K8
180 k2 = K2*XX1/K8
181 kmin2 = Kmin2/K8
182 k3 = K3/K8
183 k4 = K4*XX1/K8
184 kmin4 = Kmin4/K8
185 k5 = K5*XX1/K8
186 kmin5 = Kmin5/K8
187 k6 = K6*XX1/K8
188 kmin6 = Kmin6/K8
189 k7 = K7*XX1/K8
190 kmin7 = Kmin7/K8
191 k8 = K8/K8
192 k9 = K9/K8
193 k10 = K10/XX1
194 k11 = K11/K8
195 k12 = K12/K8
196 #
197 #------------HR--------------
198 p1 = P1*XX1/K8
199 pmin1 = Pmin1/K8
200 p2 = P2/K8
201 pmin2 = Pmin2/K8
202 p3 = P3*XX1/K8
203 pmin3 = Pmin3/K8
204 p4 = P4/K8
205 p5 = P5*XX1/K8
206 pmin5 = Pmin5/K8
207 p6 = P6*XX1/K8
208 pmin6 = Pmin6/K8
209 p7 = P7/K8
210 p8 = P8*XX1/K8
211 pmin8 = Pmin8/K8
212 p9 = P9/K8
213 p10 = P10/K8
214 p11 = P11/K8
215 p12 = P12/K8
216 #
217 #------------SSA--------------
218 q1= Q1*XX1/K8
219 qmin1 = Qmin1/K8
220 q2 = Q2*XX1/K8
221 q3 = Q3*XX1/K8
222 qmin3 = Qmin3/K8
223 q4 = Q4/K8
224 q5 = Q5*XX1/K8
225 qmin5 = Qmin5/K8
226 q6 = Q6/K8
227 #
228 #-------alt-NHEJ (MMEJ)--------
229 r1 = R1*XX1/K8
230 rmin1 = Rmin1/K8
231 r2 = R2*XX1/K8
232 r3 = R3/K8
233 r4 = R4*XX1/K8
234 rmin4 = Rmin4/K8
235 r5 = R5/K8
236 #------------------------------------
237
238 p = [ k1, k2, k3, k4, k5, k6, k7, k8, k9, k10,
239 kmin1, kmin2, kmin4, kmin5, kmin6, kmin7
240 p1, p2, p3, p4, p5, p6, p7, p8, p9, p10,
241 pmin1, pmin2, pmin3, pmin5, pmin6, pmin8
242 q1, q2, q3, q4, q5, q6,
243 qmin1, qmin3, qmin5,
244 r1, r2, r3, r4, r5,
245 rmin1, rmin4
246 ]
247
248
249 # Model definition
250 def DSBRepairPathways( Y , t , p ):
251 # SYSTEM OF DIFFERENTIAL EQUATIONS
252
253 # Concentrations of repair enzymes set to
254 # X1 # [Ku] # X2 # [DNAPKcsArt] #
255 # X4 # [PNKP] # X5 # [Pol] #
256 # X7 # [MRN/CtIP/ExoI/Dna2] #
257 # X9 # [RPA] # X10 # [Rad51/Rad51par/
258 # X11 # [DNAinc] # X12 # [Rad52] #
259 # X14 # [LigIII] # X15 # [PARP1] #
260 # X17 # [LigI]
261
262 X1= X2= X3= X4= X5= X6= X7= X8= X9= X10= X
263
264 # ----- NHEJ ----------
265 YP[0] = fNirrep - k1*Y[0]*X1 + kmin1*Y[1]
266 YP[1] = k1*Y[0]*X1 - kmin1*Y[1] - k2*Y[1]*
267 YP[2] = k2*Y[1]*X2 - k3*Y[2] - kmin2*Y[2]
268 YP[3] = k3*Y[2] - k4*(Y[3]*Y[3]) + kmin4*Y
269 YP[4] = k4*(Y[3]*Y[3]) - kmin4*Y[4] - k5*Y
270 YP[5] = kmin6*Y[6] + k5*Y[4]*X3 - kmin5*Y[
271 YP[6] = -kmin6*Y[6] - k7*Y[6]*X5 + kmin7*
272 YP[7] = k7*Y[6]*X5 - k8*Y[7] - kmin7*Y[7]
273 YP[8] = r5*Y[28] + k8*Y[7] + p12*Y[18] + p
274 YP[9] = (k9* (Y[3]+Y[4]+Y[5]+Y[6]+Y[7]) *
275
276 # ----- HR ----------
277 YP[10] = p1*Y[0]*X7 - pmin1*Y[10] - p3*Y[1
278 YP[11] = p2*X8 - pmin2*Y[11] - p3*Y[10]*Y[
279 YP[12] = p3*Y[10]*Y[11] - p4*Y[12] - pmin3
280 YP[13] = rmin1*Y[25] + p4*Y[12] - r1*X15*Y
281 YP[14] = pmin6*Y[15] + p5*Y[13]*X9 - pmin5
282 YP[15] = -p7*Y[15] - pmin6*Y[15] + p6*Y[14
283 YP[16] = p7*Y[15] - p8*Y[16]*X11 + pmin8*Y
284 YP[17] = p8*Y[16]*X11 - p9*Y[17] - pmin8*Y
285 YP[18] = p9*Y[17] - p10*Y[18] - p12*Y[18]
286 YP[19] = p10*Y[18] - p11*Y[19] # [dHJ]
287
288 # ----- SSA ----------
289 YP[20] = q1*Y[14]*X12 - qmin1*Y[20] - q2*
290 YP[21] = q2*(Y[20]*Y[20]) - q3*Y[21]*X13 +
291 YP[22] = q3*Y[21]*X13 - q4*Y[22] - qmin3*Y
292 YP[23] = q4*Y[22] - q5*Y[23]*X14 + qmin5*Y
293 YP[24] = q5*Y[23]*X14 - q6*Y[24] - qmin5*Y
294
295 # ----- MMEJ ----------
296 YP[25] = -rmin1*Y[25] - r2*Y[25]*X16 + r1*
297 YP[26] = r2*Y[25]*X16 - r3*Y[26] # [ssDNA
298 YP[27] = r3*Y[26] - r4*Y[27]*X17 + rmin4*Y
299 YP[28] = r4*Y[27]*X17 - r5*Y[28] - rmin4*Y
300
301 #-----------------------------------------
302
303 return YP
304
305 # Model execution
306 steps = od( DSBRepairPathways, Y , t , args=(p
307
308 print("\nCalculation concluded. Check the fina
309
310 # Plotting results
311 KU = []
312 DNAPKcs = []
313 RPA = []
314 Rad51 = []
315 gH2AX = []
316 gH2AXnorm = []
317 i = 0
318 for i in range(len(t)):
319 KU.append(steps[i][1])
320 DNAPKcs.append(steps[i][3]+steps[i][4]+ste
321 RPA.append(steps[i][14]+steps[i][15]+steps
322 Rad51.append(steps[i][15]+steps[i][16]+ste
323 gH2AX.append(steps[i][9])
324 gH2AXnorm.append(100*steps[i][9]/repDSB)
325
326
327 plt.figure("KU")
328 plt.plot(t, KU)
329 plt.show()
330
331
332 plt.figure("DNAPKcs")
333 plt.plot(t, DNAPKcs)
334 plt.show()
335
336
337 plt.figure("RPA")
338 plt.plot(t, RPA)
339 plt.show()
340
341
342 plt.figure("Rad51")
343 plt.plot(t, Rad51)
344 plt.show()
345
346
347 #plt.xlim(0,25)
348 plt.figure("gH2AX")
349 plt.plot(t, gH2AX)
350 plt.savefig("g-H2AX.pdf")
351
352 with open(outputFile, 'w') as ofile:
353 for m in range(len(gH2AX)):
354 ofile.write(str(t[m])+ "\t" + str(gH2A
355
356 plt.show()
357
358
359 # Normalization, and plotting with available d
360 #plt.xlim(0,25)
361 plt.figure("gH2AXnorm")
362 fig = plt.figure("gH2AXnorm")
363
364 m = max(gH2AX)
365 for k in range(len(gH2AX)):
366 gH2AXnorm[k] = 100 * gH2AX[k]/m
367
368
369 ax1 = fig.add_subplot(111)
370
371 ax1.scatter(t, gH2AXnorm, s=1, c='b', marker="
372 xD = [ 23.90410959, 10, 5,
373 yD = [ 10.55555556, 15.13888889, 26.52777778,
374 ax1.scatter(xD, yD, s=20, c='r', marker="o", l
375 xA = [ 0.5308219178082174, 2.0376712328767113
376 yA = [ 99.99999999999997, 55.13888888888887,
377 ax1.scatter(xA, yA, s=30, c='g', marker="+", l
378 plt.legend(loc='upper right');
379
380 #plt.plot(t, gH2AXnorm)
381 plt.savefig("g-H2AX_norm.pdf")
382
383 plt.show()
384
385 with open(outputFile, 'w') as ofile:
386 for m in range(len(gH2AX)):
387 ofile.write(str(t[m])+ "\t" + str(gH2A
388
389 #plt.xlim([1e1, 1e4])
390 #plt.ylim([1e-15, 3e-11])
391
392 print("\nCheck your folder for new images (pdf