| 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 molecularDNArepair.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 # 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 * b * c / Chromosome size, as defined in the mac file, but in meters. If sphere, a=b=c
19 NBP = 6405886128 # bp / Length of the DNA chain in Mbp
20 mass = 997 * 4 * 3.141592 * r3 / 3 # waterDensity * 4/3 * pi * r3 in kg
21
22 # Integration parameters
23 t0 = 0.0; # Starting time point (dimensionless)
24 t1 = 25; # Final time point (dimensionless)
25 dt = 0.001 #0.001; # Intergration time step (dimensionless)
26 numpoints = int( (t1 - t0) / dt )
27
28
29 # ONLY ADVANCED USERS AFTER THIS POINT
30 # Repair model parameters may be find at line 105-172
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.TypeClassification=="DSB+" or e.TypeClassification=="DSB++"):
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 could also be used
59
60 ffTree = f.Get("tuples/chromosome_hits")
61 for ev in ffTree:
62 acc_edep += (ev.e_chromosome_kev + ev.e_dna_kev) *1e3 # eV
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 DSB :", totalDSB, "\nRepairable DSBs :", repDSB, "\nIrrepairable DSBs :", irrDSB)
74
75 NBP_of_cell = 5000000000 # This is a constant of the model, as defined by Belov et al.
76 NBPcell = NBP/NBP_of_cell
77 totalDSBperGyCell = totalDSB/NBPcell/dose
78 print ("Dose :", dose, "\nTotal DSB/Gy/Cell :", totalDSBperGyCell)
79 S1 = repDSB/NBPcell/dose
80 S2 = irrDSB/NBPcell/dose
81 print("Repairable DSBs/Gy/cell :", S1, "\nIrrepairable DSBs/Gy/cell :", S2)
82
83
84 # Model parameters
85 fDz = 1 # Dose is set to 1 Gy, because falpha is normalized to DSB/Gy/Cell
86 falpha = totalDSBperGyCell #27.5 #totalDSB # total DSB
87 fNirrep = irrDSB/totalDSB # irrepDSB # racirreparable = (G4double)fIrreparableDSBs/(G4double)fTotDSBs;
88 if (fNirrep == 0):
89 fNirrep = 0.01 # For not printing error
90 print("Ratio of irrepairable to repairable damage :", fNirrep)
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 the ODE solver.
101 # I use a large number of points, only because I want to make
102 # a plot of the solution that looks nice.
103 t = [(t1-t0) * float(i) / (numpoints - 1) for i in range(numpoints)]
104
105
106 # DIMENSIONAL REACTION RATES DEFINITON (check Belov et al.) (Belov OV, et al. J Theor Biol. 2015 Feb 7;366:115-30. doi: 10.1016/j.jtbi.2014.09.024)
107 #
108 #------------NHEJ--------------
109 K1 = 11.052 #other:10.02 #11.052 # M-1*h-1 # change parametrs based on the article of Belov et al.
110 Kmin1 = 6.6*1e-01 #6.59999*1e-04 # h-1 # other 6.6*1e-01 #
111 K2 = 5.82*1e+05 #18.8305*(1.08517-np.exp(-21.418/pow(fDz,1.822))) # M-1*h-1 # other:5.82*1e+05 #
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 needs to be adjusted.
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-1
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*h-1
153 Qmin1 = 1.71*1e-04 # h-1
154 Q2 = 3*1e+04 #4.8052*1e+04 # M-1*h-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, k11, k12,
239 kmin1, kmin2, kmin4, kmin5, kmin6, kmin7,
240 p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12,
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 be constant
254 # X1 # [Ku] # X2 # [DNAPKcsArt] # X3 # [LigIV/XRCC4/XLF]
255 # X4 # [PNKP] # X5 # [Pol] # X6 # [H2AX]
256 # X7 # [MRN/CtIP/ExoI/Dna2] # X8 # [ATM]
257 # X9 # [RPA] # X10 # [Rad51/Rad51par/BRCA2]
258 # X11 # [DNAinc] # X12 # [Rad52] # X13 # [ERCC1/XPF]
259 # X14 # [LigIII] # X15 # [PARP1] # X16 # [Pol]
260 # X17 # [LigI]
261
262 X1= X2= X3= X4= X5= X6= X7= X8= X9= X10= X11= X12= X13= X14= X15= X16= X17= 400000.
263
264 # ----- NHEJ ----------
265 YP[0] = fNirrep - k1*Y[0]*X1 + kmin1*Y[1] - p1*Y[0]*X1 + pmin1*Y[10] # [DSB]
266 YP[1] = k1*Y[0]*X1 - kmin1*Y[1] - k2*Y[1]*X2 + kmin2*Y[2] # [DBS * Ku]
267 YP[2] = k2*Y[1]*X2 - k3*Y[2] - kmin2*Y[2] # [DSB * DNA-PK/Art]
268 YP[3] = k3*Y[2] - k4*(Y[3]*Y[3]) + kmin4*Y[4] # [DSB * DNA-PK/ArtP]
269 YP[4] = k4*(Y[3]*Y[3]) - kmin4*Y[4] - k5*Y[4]*X3 + kmin5*Y[5] # [Bridge]
270 YP[5] = kmin6*Y[6] + k5*Y[4]*X3 - kmin5*Y[5] - k6*Y[5]*X4 # [Bridge * LigIV/XRCC4/XLF]
271 YP[6] = -kmin6*Y[6] - k7*Y[6]*X5 + kmin7*Y[7] + k6*Y[5]*X4 # [Bridge * LigIV/XRCC4/XLF * PNKP]
272 YP[7] = k7*Y[6]*X5 - k8*Y[7] - kmin7*Y[7] # [Bridge * LigIV/XRCC4/XLF * PNKP * Pol]
273 YP[8] = r5*Y[28] + k8*Y[7] + p12*Y[18] + p11*Y[19] + q6*Y[24] # [dsDNA]
274 YP[9] = (k9* (Y[3]+Y[4]+Y[5]+Y[6]+Y[7]) * X6)/ (k10+ Y[3]+ Y[4]+ Y[5]+ Y[6]+ Y[7]) - k11*Y[8]- k12*Y[9] # [gH2AX foci]
275
276 # ----- HR ----------
277 YP[10] = p1*Y[0]*X7 - pmin1*Y[10] - p3*Y[10]*Y[11] + pmin3*Y[12] # [MRN/CtIP/ExoI/Dna2]
278 YP[11] = p2*X8 - pmin2*Y[11] - p3*Y[10]*Y[11] + p4*Y[12] + pmin3*Y[12] # [ATMP]
279 YP[12] = p3*Y[10]*Y[11] - p4*Y[12] - pmin3*Y[12] # [DSB * MRN/CtIP/ExoI/Dna2 * ATMP]
280 YP[13] = rmin1*Y[25] + p4*Y[12] - r1*X15*Y[13] - p5*Y[13]*X9 + pmin5*Y[14] # [ssDNA]
281 YP[14] = pmin6*Y[15] + p5*Y[13]*X9 - pmin5*Y[14] - p6*Y[14]*X10 - q1*Y[14]*X12 + qmin1*Y[20] # [ssDNA * RPA]
282 YP[15] = -p7*Y[15] - pmin6*Y[15] + p6*Y[14]*X10 # [ssDNA * RPA * Rad51/Rad51par/BRCA2]
283 YP[16] = p7*Y[15] - p8*Y[16]*X11 + pmin8*Y[17] # [Rad51 filament]
284 YP[17] = p8*Y[16]*X11 - p9*Y[17] - pmin8*Y[17] # [Rad51 filament * DNAinc]
285 YP[18] = p9*Y[17] - p10*Y[18] - p12*Y[18] # [D-loop]
286 YP[19] = p10*Y[18] - p11*Y[19] # [dHJ]
287
288 # ----- SSA ----------
289 YP[20] = q1*Y[14]*X12 - qmin1*Y[20] - q2*(Y[20]*Y[20]) # [ssDNA * RPA * Rad52]
290 YP[21] = q2*(Y[20]*Y[20]) - q3*Y[21]*X13 + qmin3*Y[22] # [Flap]
291 YP[22] = q3*Y[21]*X13 - q4*Y[22] - qmin3*Y[22] # [Flap * ERCC1/XPF]
292 YP[23] = q4*Y[22] - q5*Y[23]*X14 + qmin5*Y[24] # [dsDNAnicks]
293 YP[24] = q5*Y[23]*X14 - q6*Y[24] - qmin5*Y[24] # [dsDNAnicks * LigIII]
294
295 # ----- MMEJ ----------
296 YP[25] = -rmin1*Y[25] - r2*Y[25]*X16 + r1*X15*Y[13] # [ssDNA * PARP1]
297 YP[26] = r2*Y[25]*X16 - r3*Y[26] # [ssDNA * Pol]
298 YP[27] = r3*Y[26] - r4*Y[27]*X17 + rmin4*Y[28] # [MicroHomol]
299 YP[28] = r4*Y[27]*X17 - r5*Y[28] - rmin4*Y[28] # [MicroHomol * LigI]
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 final images and close them to end the program.!")
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]+steps[i][5]+steps[i][6]+steps[i][7])
321 RPA.append(steps[i][14]+steps[i][15]+steps[i][20])
322 Rad51.append(steps[i][15]+steps[i][16]+steps[i][17])
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(gH2AX[m])+"\n")
355
356 plt.show()
357
358
359 # Normalization, and plotting with available data, as published in Chatzipapas et al., arxiv https://doi.org/10.48550/arXiv.2210.01564
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="s", label='gH2AXnorm')
372 xD = [ 23.90410959, 10, 5, 1.47260274, 0.71917808, 0.239726027, 0.136986301, 0.102739726 ]
373 yD = [ 10.55555556, 15.13888889, 26.52777778, 80, 100, 60, 20, 8.333333333 ]
374 ax1.scatter(xD, yD, s=20, c='r', marker="o", label='gH2AXnorm_dous')
375 xA = [ 0.5308219178082174, 2.0376712328767113, 4.006849315068493, 11.952054794520548, 23.869863013698634 ]
376 yA = [ 99.99999999999997, 55.13888888888887, 30.27777777777777, 12.777777777777771, 1.6666666666666572 ]
377 ax1.scatter(xA, yA, s=30, c='g', marker="+", label='gH2AXnorm_Asaithamby')
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(gH2AXnorm[m])+"\n")
388
389 #plt.xlim([1e1, 1e4])
390 #plt.ylim([1e-15, 3e-11])
391
392 print("\nCheck your folder for new images (pdf). Thank you for using molecularDNArepair.!\n")