moleculardna/repair_survival_models/molecularDNArepair.py

0001 #!/usr/bin/env python
0002 # coding: utf-8
0003
0004 # To run, change the links and definitions in this script and
0005 # type in the terminal containing the file
0006 # $python3 molecularDNArepair.py
0007
0008 # Authors:
0009 # K. Chatzipapas (*), D. Sakata, H. N. Tran, M. Dordevic, S. Incerti et al.
0010 # (*) contact: k.chatzipapas@yahoo.com
0011
0012 # Define filenames
0013 outputFile = "molecularDNArepair.txt"
0014 iRootFile  = "../molecular-dna.root"
0015 print("\nInput file:",iRootFile,"\n")
0016
0017 # Define cell parameters
0018 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
0019 NBP = 6405886128 # bp / Length of the DNA chain in Mbp
0020 mass = 997 * 4 * 3.141592 * r3 / 3     # waterDensity * 4/3 * pi * r3 in kg
0021
0022 # Integration parameters
0023 t0 = 0.0;   # Starting time point (dimensionless)
0024 t1 = 25;  # Final time point (dimensionless)
0025 dt = 0.001 #0.001; # Intergration time step (dimensionless)
0026 numpoints = int( (t1 - t0) / dt )
0027
0028
0029 # ONLY ADVANCED USERS AFTER THIS POINT
0030 # Repair model parameters may be find at line 105-172
0031 #####################################################################################################################
0032 import ROOT as root
0033 import numpy as np
0034 import matplotlib.pyplot as plt
0035 from scipy.integrate import odeint as od
0036
0037 # Initialization
0038 DSBBPID = []
0039 i = 0
0040 repDSB = 0
0041 irrDSB = 0
0042 acc_edep = 0
0043 dose = 0
0044 eVtoJ = 1.60218e-19
0045
0046
0047 f = root.TFile(iRootFile)
0048
0049 fTree = f.Get("tuples/damage")
0050 for e in fTree:
0051     if (e.TypeClassification=="DSB" or e.TypeClassification=="DSB+" or e.TypeClassification=="DSB++"):
0052         DSBBPID.append((e.Event, e.BasePair))
0053         i += 1
0054
0055 gTree = f.Get("tuples/classification")
0056 for e in gTree:
0057     repDSB += e.DSB
0058     irrDSB += e.DSBp + e.DSBpp  # 2*eDSBpp   could also be used
0059
0060 ffTree = f.Get("tuples/chromosome_hits")
0061 for ev in ffTree:
0062     acc_edep += (ev.e_chromosome_kev + ev.e_dna_kev) *1e3  # eV
0063
0064
0065
0066 # Calculate the total number of DSBs
0067 totalDSB = repDSB + irrDSB
0068
0069 # Calculate the absorbed dose
0070 dose = acc_edep * eVtoJ / mass    # Dose in Gy
0071
0072 #print ("r3=", r3)
0073 #print("Dose :", dose, "\ni :", i, "\nTotal DSB :", totalDSB, "\nRepairable DSBs :", repDSB, "\nIrrepairable DSBs :", irrDSB)
0074
0075 NBP_of_cell = 5000000000   # This is a constant of the model, as defined by Belov et al.
0076 NBPcell = NBP/NBP_of_cell
0077 totalDSBperGyCell = totalDSB/NBPcell/dose
0078 print ("Dose :", dose, "\nTotal DSB/Gy/Cell :", totalDSBperGyCell)
0079 S1 = repDSB/NBPcell/dose
0080 S2 = irrDSB/NBPcell/dose
0081 print("Repairable DSBs/Gy/cell :", S1, "\nIrrepairable DSBs/Gy/cell :", S2)
0082
0083
0084 # Model parameters
0085 fDz     = 1 # Dose is set to 1 Gy, because falpha is normalized to DSB/Gy/Cell
0086 falpha  = totalDSBperGyCell #27.5 #totalDSB # total DSB
0087 fNirrep = irrDSB/totalDSB # irrepDSB # racirreparable = (G4double)fIrreparableDSBs/(G4double)fTotDSBs;
0088 if (fNirrep == 0):
0089     fNirrep = 0.01   # For not printing error
0090 print("Ratio of irrepairable to repairable damage :", fNirrep)
0091
0092 # Initial conditions for ODE
0093 NbEquat = 29
0094 Y  = [0] * NbEquat
0095 #print (Y[1])
0096 Y[0] = falpha * fDz
0097 YP = [0] * NbEquat
0098 print("\nY[0] used in the model :", Y[0])
0099
0100 # Create the time samples for the output of the ODE solver.
0101 # I use a large number of points, only because I want to make
0102 # a plot of the solution that looks nice.
0103 t = [(t1-t0) * float(i) / (numpoints - 1) for i in range(numpoints)]
0104
0105
0106 #  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)
0107 #
0108 #------------NHEJ--------------
0109 K1 = 11.052 #other:10.02         #11.052             # M-1*h-1  # change parametrs based on the article of Belov et al.
0110 Kmin1 = 6.6*1e-01       #6.59999*1e-04   # h-1   # other 6.6*1e-01       #
0111 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    #
0112 #print (K2)
0113 Kmin2 = 5.26*1e-01      #h-1
0114 K3 = 1.86               # h-1
0115 K4 = 1.38*1e+06         # M-1*h-1
0116 Kmin4 = 3.86*1e-04      # h-1
0117 K5 = 15.24              # M-1*h-1
0118 Kmin5 = 8.28            # h-1
0119 K6 = 18.06              # M-1*h-1
0120 Kmin6 = 1.33            # h-1
0121 K7 = 2.73*1e+05         # M-1*h-1
0122 Kmin7 = 3.2             # h-1
0123 #K8 is a normalization factor and sometimes it needs to be adjusted.
0124 normK8 = 0.35
0125 K8 = normK8*5.52*1e-01  # h-1
0126 K9 = 1.66*1e-01         # h-1
0127 K10 = (1.93*1e-07)/fNirrep # M
0128 K11 = 7.50*1e-02        # h-1
0129 K12 = 11.1              # h-1
0130 #
0131 #------------HR--------------
0132 P1 = 1.75*1e+03         # M-1*h-1
0133 Pmin1 = 1.33*1e-04      # h-1
0134 P2 = 0.39192            # h-1
0135 Pmin2 = 2.7605512*1e+02 # h-1
0136 P3 = 1.37*1e+04         # M-1*h-1
0137 Pmin3 = 2.34            # h-1
0138 P4 = 5.52*1e-2     #3.588*1e-02        # h-1
0139 P5 = 1.20*1e+05         # M-1*h-1
0140 Pmin5 = 8.82*1e-05      # h-1
0141 P6 = 1.87*1e+05    #1.54368*1e+06      # M-1*h-1
0142 Pmin6 = 1.55*1e-03      # h-1
0143 P7 = 21.36         #1.4904             # h-1
0144 P8 = 1.20*1e+04         # M-1*h-1
0145 Pmin8 = 2.49*1e-04      # h-1
0146 P9 = 4.88*1e-1     #1.104              # h-1
0147 P10 = 7.20*1e-03        # h-1
0148 P11 = 6.06*1e-04        # h-1
0149 P12 = 2.76*1e-01        # h-1
0150 #
0151 #------------SSA--------------
0152 Q1 = 7.8*1e+03      #1.9941*1e+05       # M-1*h-1
0153 Qmin1 = 1.71*1e-04      # h-1
0154 Q2 = 3*1e+04        #4.8052*1e+04       # M-1*h-1
0155 Q3 = 6*1e+03            # M-1*h-1
0156 Qmin3 = 6.06*1e-04      # h-1
0157 Q4 = 1.66*1e-06     #1.62*1e-03         # h-1
0158 Q5 = 8.40*1e+04         # M-1*h-1
0159 Qmin5 = 4.75*1e-04      # h-1
0160 Q6 = 11.58              # h-1
0161 #
0162 #-------alt-NHEJ (MMEJ)--------
0163 R1 = 2.39*1e+03         # M-1*h-1
0164 Rmin1 = 12.63           # h-1
0165 R2 = 4.07*1e+04         # M-1*h-1
0166 R3 = 9.82               # h-1
0167 R4 = 1.47*1e+05         # M-1*h-1
0168 Rmin4 = 2.72            # h-1
0169 R5 = 1.65*1e-01         # h-1
0170 #
0171 # Scalling rate XX1
0172 XX1 = 9.19*1e-07 # M
0173 XX3 = 2.3 *1e-12 # M
0174 #
0175 # DIMENSIONLESS REACTION RATES
0176 #
0177 #------------NHEJ--------------
0178 k1 = K1*XX1/K8
0179 kmin1 = Kmin1/K8
0180 k2 = K2*XX1/K8
0181 kmin2 = Kmin2/K8
0182 k3 = K3/K8
0183 k4 = K4*XX1/K8
0184 kmin4 = Kmin4/K8
0185 k5 = K5*XX1/K8
0186 kmin5 = Kmin5/K8
0187 k6 = K6*XX1/K8
0188 kmin6 = Kmin6/K8
0189 k7 = K7*XX1/K8
0190 kmin7 = Kmin7/K8
0191 k8 = K8/K8
0192 k9 = K9/K8
0193 k10 = K10/XX1
0194 k11 = K11/K8
0195 k12 = K12/K8
0196 #
0197 #------------HR--------------
0198 p1 = P1*XX1/K8
0199 pmin1 = Pmin1/K8
0200 p2 = P2/K8
0201 pmin2 = Pmin2/K8
0202 p3 = P3*XX1/K8
0203 pmin3 = Pmin3/K8
0204 p4 = P4/K8
0205 p5 = P5*XX1/K8
0206 pmin5 = Pmin5/K8
0207 p6 = P6*XX1/K8
0208 pmin6 = Pmin6/K8
0209 p7 = P7/K8
0210 p8 = P8*XX1/K8
0211 pmin8 = Pmin8/K8
0212 p9 = P9/K8
0213 p10 = P10/K8
0214 p11 = P11/K8
0215 p12 = P12/K8
0216 #
0217 #------------SSA--------------
0218 q1= Q1*XX1/K8
0219 qmin1 = Qmin1/K8
0220 q2 = Q2*XX1/K8
0221 q3 = Q3*XX1/K8
0222 qmin3 = Qmin3/K8
0223 q4 = Q4/K8
0224 q5 = Q5*XX1/K8
0225 qmin5 = Qmin5/K8
0226 q6 = Q6/K8
0227 #
0228 #-------alt-NHEJ (MMEJ)--------
0229 r1 = R1*XX1/K8
0230 rmin1 = Rmin1/K8
0231 r2 = R2*XX1/K8
0232 r3 = R3/K8
0233 r4 = R4*XX1/K8
0234 rmin4 = Rmin4/K8
0235 r5 = R5/K8
0236 #------------------------------------
0237
0238 p = [ k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12,
0239       kmin1, kmin2, kmin4, kmin5, kmin6, kmin7,
0240       p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12,
0241       pmin1, pmin2, pmin3, pmin5, pmin6, pmin8,
0242       q1, q2, q3, q4, q5, q6,
0243       qmin1, qmin3, qmin5,
0244       r1, r2, r3, r4, r5,
0245       rmin1, rmin4
0246       ]
0247
0248
0249 # Model definition
0250 def DSBRepairPathways( Y , t , p ):
0251     # SYSTEM OF DIFFERENTIAL EQUATIONS
0252
0253     # Concentrations of repair enzymes set to be constant
0254     # X1  # [Ku]      # X2  # [DNAPKcsArt]   # X3  # [LigIV/XRCC4/XLF]
0255     # X4  # [PNKP]    # X5  # [Pol]          # X6  # [H2AX]
0256     # X7  # [MRN/CtIP/ExoI/Dna2]             # X8  # [ATM]
0257     # X9  # [RPA]     # X10 # [Rad51/Rad51par/BRCA2]
0258     # X11 # [DNAinc]  # X12 # [Rad52]        # X13 # [ERCC1/XPF]
0259     # X14 # [LigIII]  # X15 # [PARP1]        # X16 # [Pol]
0260     # X17 # [LigI]
0261
0262     X1= X2= X3= X4= X5= X6= X7= X8= X9= X10= X11= X12= X13= X14= X15= X16= X17=  400000.
0263
0264     # ----- NHEJ ----------
0265     YP[0] = fNirrep - k1*Y[0]*X1 + kmin1*Y[1] - p1*Y[0]*X1 + pmin1*Y[10] # [DSB]
0266     YP[1] = k1*Y[0]*X1 - kmin1*Y[1] - k2*Y[1]*X2 + kmin2*Y[2] # [DBS * Ku]
0267     YP[2] = k2*Y[1]*X2 - k3*Y[2] - kmin2*Y[2] # [DSB * DNA-PK/Art]
0268     YP[3] = k3*Y[2] - k4*(Y[3]*Y[3]) + kmin4*Y[4] # [DSB * DNA-PK/ArtP]
0269     YP[4] = k4*(Y[3]*Y[3]) - kmin4*Y[4] - k5*Y[4]*X3 + kmin5*Y[5] # [Bridge]
0270     YP[5] = kmin6*Y[6] + k5*Y[4]*X3 - kmin5*Y[5] - k6*Y[5]*X4 # [Bridge * LigIV/XRCC4/XLF]
0271     YP[6] = -kmin6*Y[6] - k7*Y[6]*X5 +  kmin7*Y[7] + k6*Y[5]*X4 # [Bridge * LigIV/XRCC4/XLF * PNKP]
0272     YP[7] = k7*Y[6]*X5 - k8*Y[7] - kmin7*Y[7] # [Bridge * LigIV/XRCC4/XLF * PNKP * Pol]
0273     YP[8] = r5*Y[28] + k8*Y[7] + p12*Y[18] + p11*Y[19] + q6*Y[24] # [dsDNA]
0274     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]
0275
0276     # ----- HR ----------
0277     YP[10] = p1*Y[0]*X7 - pmin1*Y[10] - p3*Y[10]*Y[11] + pmin3*Y[12]     # [MRN/CtIP/ExoI/Dna2]
0278     YP[11] = p2*X8 - pmin2*Y[11] - p3*Y[10]*Y[11] + p4*Y[12] + pmin3*Y[12]     # [ATMP]
0279     YP[12] = p3*Y[10]*Y[11] - p4*Y[12] - pmin3*Y[12]     # [DSB * MRN/CtIP/ExoI/Dna2 * ATMP]
0280     YP[13] = rmin1*Y[25] + p4*Y[12] - r1*X15*Y[13] - p5*Y[13]*X9 + pmin5*Y[14]     # [ssDNA]
0281     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]
0282     YP[15] = -p7*Y[15] - pmin6*Y[15] + p6*Y[14]*X10    # [ssDNA * RPA * Rad51/Rad51par/BRCA2]
0283     YP[16] = p7*Y[15] - p8*Y[16]*X11 + pmin8*Y[17] # [Rad51 filament]
0284     YP[17] = p8*Y[16]*X11 - p9*Y[17] - pmin8*Y[17] # [Rad51 filament * DNAinc]
0285     YP[18] = p9*Y[17] - p10*Y[18] - p12*Y[18] # [D-loop]
0286     YP[19] = p10*Y[18] - p11*Y[19] # [dHJ]
0287
0288     # ----- SSA ----------
0289     YP[20] = q1*Y[14]*X12 - qmin1*Y[20] -  q2*(Y[20]*Y[20])    # [ssDNA * RPA * Rad52]
0290     YP[21] = q2*(Y[20]*Y[20]) - q3*Y[21]*X13 + qmin3*Y[22] # [Flap]
0291     YP[22] = q3*Y[21]*X13 - q4*Y[22] - qmin3*Y[22] # [Flap * ERCC1/XPF]
0292     YP[23] = q4*Y[22] - q5*Y[23]*X14 + qmin5*Y[24] # [dsDNAnicks]
0293     YP[24] = q5*Y[23]*X14 - q6*Y[24] - qmin5*Y[24] # [dsDNAnicks * LigIII]
0294
0295     # ----- MMEJ ----------
0296     YP[25] = -rmin1*Y[25] - r2*Y[25]*X16 + r1*X15*Y[13] # [ssDNA * PARP1]
0297     YP[26] = r2*Y[25]*X16 - r3*Y[26] # [ssDNA * Pol]
0298     YP[27] = r3*Y[26] - r4*Y[27]*X17 + rmin4*Y[28] # [MicroHomol]
0299     YP[28] = r4*Y[27]*X17 - r5*Y[28] - rmin4*Y[28] # [MicroHomol * LigI]
0300
0301     #---------------------------------------------------
0302
0303     return YP
0304
0305 # Model execution
0306 steps = od( DSBRepairPathways, Y , t , args=(p,) )
0307
0308 print("\nCalculation concluded. Check the final images and close them to end the program.!")
0309
0310 # Plotting results
0311 KU = []
0312 DNAPKcs = []
0313 RPA = []
0314 Rad51 = []
0315 gH2AX = []
0316 gH2AXnorm = []
0317 i = 0
0318 for i in range(len(t)):
0319     KU.append(steps[i][1])
0320     DNAPKcs.append(steps[i][3]+steps[i][4]+steps[i][5]+steps[i][6]+steps[i][7])
0321     RPA.append(steps[i][14]+steps[i][15]+steps[i][20])
0322     Rad51.append(steps[i][15]+steps[i][16]+steps[i][17])
0323     gH2AX.append(steps[i][9])
0324     gH2AXnorm.append(100*steps[i][9]/repDSB)
0325
0326
0327 plt.figure("KU")
0328 plt.plot(t, KU)
0329 plt.show()
0330
0331
0332 plt.figure("DNAPKcs")
0333 plt.plot(t, DNAPKcs)
0334 plt.show()
0335
0336
0337 plt.figure("RPA")
0338 plt.plot(t, RPA)
0339 plt.show()
0340
0341
0342 plt.figure("Rad51")
0343 plt.plot(t, Rad51)
0344 plt.show()
0345
0346
0347 #plt.xlim(0,25)
0348 plt.figure("gH2AX")
0349 plt.plot(t, gH2AX)
0350 plt.savefig("g-H2AX.pdf")
0351
0352 with open(outputFile, 'w') as ofile:
0353     for m in range(len(gH2AX)):
0354         ofile.write(str(t[m])+ "\t" + str(gH2AX[m])+"\n")
0355
0356 plt.show()
0357
0358
0359 # Normalization, and plotting with available data, as published in Chatzipapas et al., arxiv https://doi.org/10.48550/arXiv.2210.01564
0360 #plt.xlim(0,25)
0361 plt.figure("gH2AXnorm")
0362 fig = plt.figure("gH2AXnorm")
0363
0364 m = max(gH2AX)
0365 for k in range(len(gH2AX)):
0366     gH2AXnorm[k] = 100 * gH2AX[k]/m
0367
0368
0369 ax1 = fig.add_subplot(111)
0370
0371 ax1.scatter(t, gH2AXnorm, s=1, c='b', marker="s", label='gH2AXnorm')
0372 xD = [ 23.90410959, 10,           5,           1.47260274, 0.71917808, 0.239726027, 0.136986301, 0.102739726 ]
0373 yD = [ 10.55555556, 15.13888889, 26.52777778, 80,        100,          60,         20,           8.333333333 ]
0374 ax1.scatter(xD, yD, s=20, c='r', marker="o", label='gH2AXnorm_dous')
0375 xA = [  0.5308219178082174, 2.0376712328767113, 4.006849315068493, 11.952054794520548, 23.869863013698634 ]
0376 yA = [ 99.99999999999997,  55.13888888888887,  30.27777777777777,  12.777777777777771,  1.6666666666666572 ]
0377 ax1.scatter(xA, yA, s=30, c='g', marker="+", label='gH2AXnorm_Asaithamby')
0378 plt.legend(loc='upper right');
0379
0380 #plt.plot(t, gH2AXnorm)
0381 plt.savefig("g-H2AX_norm.pdf")
0382
0383 plt.show()
0384
0385 with open(outputFile, 'w') as ofile:
0386     for m in range(len(gH2AX)):
0387         ofile.write(str(t[m])+ "\t" + str(gH2AXnorm[m])+"\n")
0388
0389 #plt.xlim([1e1, 1e4])
0390 #plt.ylim([1e-15, 3e-11])
0391
0392 print("\nCheck your folder for new images (pdf). Thank you for using molecularDNArepair.!\n")