moleculardna/repair_survival_models/molecularDNAsurvival.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 molecularDNAsurvival.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 # Authors' note : This is a very early version of this calculation technique.
0013 #                 Needs to optimized for each different application. An updated
0014 #                 version will be released in the future.
0015
0016
0017 # Define filenames
0018 outputFile = "molecularDNAsurvival.txt"
0019 iRootFile  = "../molecular-dna.root"
0020 print("\nInput file:",iRootFile,"\n")
0021
0022 # Define cell parameters
0023 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
0024 NBP = 6405886128 # bp / Length of the DNA chain in Mbp
0025 mass = 997 * 4 * 3.141592 * r3 / 3     # waterDensity * 4/3 * pi * r3 in kg
0026
0027 # Integration parameters
0028 t0 = 0.0   # Starting time point (dimensionless)
0029 t1 = 336   # Final time point (dimensionless)
0030 dt = 0.001 # Intergration time step (dimensionless)
0031 numpoints = int( (t1 - t0) / dt )
0032
0033 # Survival model Parameters
0034 # THESE PARAMETERS HAVE TO BE DEFINED FOR PROPER RESULTS/ EXPERIMENTAL PARAMETERS,
0035 # AS DESCRIBED IN (Stewart RD. Radiat Res. 2001 https://pubmed.ncbi.nlm.nih.gov/11554848/ )
0036 maxDose = 9
0037 DR1     = 720      #Gy/h SF-dose
0038 bp      = 1        #Gbp
0039 #DR2   = 600.      #Gy/h SF-time
0040 #DR3   = 12000.    #Gy/h FAR
0041
0042 # usual value of gamma 0.25 (changes the end value of the curve)
0043 gamma = 0.19  #other published values: 0.189328
0044 # Lamb1 the end of the curve Lamb2 is connected with Eta
0045 Lamb1 = 0.10  #other published values: 0.0125959   # 0.671  h-1    or 2.77 h-1 (15 min repair half-time)
0046 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)
0047 #
0048 Beta1 = 0     #other published values: 0.0193207   # 0.00152
0049 Beta2 = 0
0050 # Defines the curvature as well as the end value of the curve
0051 Eta   = 0.00005  #other published values: #7.50595e-06 # 1.18e-04 h-1  or 0.00042 h-1
0052 # name of the cell (used in the output)
0053 cell  = "test"
0054
0055 # ONLY ADVANCED USERS AFTER THIS POINT
0056 #####################################################################################################################
0057 import ROOT as root
0058 import numpy as np
0059 import matplotlib.pyplot as plt
0060 from scipy.integrate import ode
0061 import threading
0062
0063 # Definition of the survival model
0064 def SF_system( T_, DR_, Y_, Sig1_, Sig2_, gamma_, Lamb1_, Lamb2_, Beta1_, Beta2_, Eta_ ):
0065
0066     global T
0067     global DR
0068     global Y
0069     global Sig1
0070     global Sig2
0071     global gamma
0072     global Lamb1
0073     global Lamb2
0074     global Beta1
0075     global Beta2
0076     global Eta
0077
0078     T  = T_
0079     DR = DR_
0080     Y  = Y_
0081     Sig1  = Sig1_
0082     Sig2  = Sig2_
0083     gamma = gamma_
0084     Lamb1 = Lamb1_
0085     Lamb2 = Lamb2_
0086     Beta1 = Beta1_
0087     Beta2 = Beta2_
0088     Eta   = Eta_
0089
0090
0091     return T, DR, Y, Sig1, Sig2, gamma, Lamb1, Lamb2, Beta1, Beta2, Eta
0092
0093 def SF_function( t, x, p ):
0094     if(t<=T):
0095         #print("t<=T")
0096         dx[0] = DR*Y*Sig1 - Lamb1*x[0] - Eta*x[0]*(x[0]+x[1])
0097         dx[1] = DR*Y*Sig2 - Lamb2*x[1] - Eta*x[1]*(x[0]+x[1])
0098         dx[2] = Beta1*Lamb1*x[0] + Beta2*Lamb2*x[1]+gamma*Eta*(x[0]+x[1])*(x[0]+x[1])
0099         dx[3] = (1.-Beta1)*Lamb1*x[0] + (1.-Beta2)*Lamb2*x[1]+(1.-gamma)*Eta*(x[0]+x[1])*(x[0]+x[1])
0100     else:
0101         dx[0] =  - Lamb1*x[0] - Eta*x[0]*(x[0]+x[1])
0102         dx[1] =  - Lamb2*x[1] - Eta*x[1]*(x[0]+x[1])
0103         dx[2] = Beta1*Lamb1*x[0] + Beta2*Lamb2*x[1]+gamma*Eta*(x[0]+x[1])*(x[0]+x[1])
0104         dx[3] = (1.-Beta1)*Lamb1*x[0] + (1.-Beta2)*Lamb2*x[1]+(1.-gamma)*Eta*(x[0]+x[1])*(x[0]+x[1])
0105     return dx
0106
0107
0108 DSBBPID = []
0109 repDSB = 0
0110 irrDSB = 0
0111 acc_edep = 0
0112 dose = 0
0113 i = 0
0114 eVtoJ = 1.60218e-19
0115
0116
0117 # Analyze root files
0118 f = root.TFile(iRootFile)
0119
0120 fTree = f.Get("tuples/damage")
0121 for e in fTree:
0122     if (e.TypeClassification=="DSB" or e.TypeClassification=="DSB+" or e.TypeClassification=="DSB++"):
0123         DSBBPID.append((e.Event, e.BasePair))
0124         i += 1
0125
0126 gTree = f.Get("tuples/classification")
0127 for e in gTree:
0128     repDSB += e.DSB
0129     irrDSB += e.DSBp + 2*e.DSBpp
0130
0131 ffTree = f.Get("tuples/chromosome_hits")
0132 for ev in ffTree:
0133     acc_edep += (ev.e_chromosome_kev + ev.e_dna_kev) *1e3  # eV
0134
0135
0136 # Calculate the total number of DSBs
0137 totalDSB = repDSB + irrDSB
0138 # Calculate the absorbed dose
0139 dose = acc_edep * eVtoJ / mass    # Dose in Gy
0140
0141 #print ("r3=", r3)
0142 #print("Dose :", dose, "\nTotal DSB/Gy :", totalDSB/dose, "\nRepairable DSBs/Gy :", repDSB/dose, "\nIrrepairable DSBs/Gy :", irrDSB/dose)
0143
0144 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/ )
0145 NBPcell = NBP/NBP_of_cell
0146 totalDSBperGyCell = totalDSB/NBPcell/dose
0147 #totalDSBperGyCell = totalDSB/dose
0148 print ("Dose :", dose, "\nTotal DSB/Gy/Cell :", totalDSBperGyCell)
0149
0150 kFactor = 1  # This factor is used to normalize irrepairable damage to lower values than 1.
0151 S1 = repDSB/NBPcell/dose/kFactor
0152 S2 = irrDSB/NBPcell/dose/kFactor
0153 print("Repairable DSBs/Gy/cell :", S1, "\nIrrepairable DSBs/Gy/cell :", S2, "\n\n")
0154
0155
0156 # In[ ]:
0157
0158
0159 T  = 0
0160 DR = 0
0161 Y  = 0
0162 kozi = False
0163 sum = [0]*(maxDose+1)
0164 outputfile = [""] * (maxDose+1)
0165
0166 # Create the time samples for the output of the ODE solver.
0167 # I use a large number of points, only because I want to make
0168 # a plot of the solution that looks nice.
0169 #t = [(t1-t0) * float(i) / (numpoints - 1) for i in range(numpoints)]
0170 t = np.linspace(t0, t1, numpoints)
0171
0172 dx = np.zeros(4)
0173 x  = np.zeros(4)
0174 dsbs = 0.0
0175 Sig1 = S1
0176 Sig2 = S2
0177
0178 filename =  "txt/molecularDNAsurvival_"+cell+".dat"
0179 outputfile1 = open(filename,"w")
0180 outputfile1.write("Dose \t Survival Fraction\n")
0181
0182 for j in range(maxDose+1):
0183     #j=1
0184     sol = []
0185     time = []
0186     dose = j
0187     expEndTime = dose/DR1 #h
0188
0189     #////double StopTime  = (expEndTime+ n/min(Lamb1,Lamb2));
0190
0191     arguments = SF_system( expEndTime, DR1,bp,Sig1,Sig2,gamma,Lamb1,Lamb2,Beta1,Beta2,Eta )
0192     p = [ T, DR, Y, Sig1, Sig2, gamma, Lamb1, Lamb2, Beta1, Beta2, Eta ]
0193
0194     Stepper = "dopri5"  # dopri5,vode,dop853
0195
0196     if (j==0):
0197         solver = ode(SF_function).set_integrator(Stepper, nsteps=1)
0198         print("#*#*#*#*#*# Ignore this warning #*#*#*#*#*#")
0199     else:
0200         solver = ode(SF_function).set_integrator(Stepper, nsteps=numpoints)
0201     solver.set_initial_value( x, t0 ).set_f_params(p)
0202     #print(solver.t,solver.y)
0203
0204     # Repeatedly call the `integrate` method to advance the
0205     # solution to time t+dt, and save the solution in solver.
0206     while solver.t < t1: # and solver.successful():
0207         solver.integrate(solver.t+dt)
0208
0209     sum[j] = solver.y[2]
0210     print("Survival probability of cell ("+cell+") after",j,"Gy is", np.exp(-sum[j]) )
0211     outputfile1.write(str(dose) +"\t"+ str(np.exp(-(solver.y[2]))) +"\n")
0212
0213 c = np.linspace(0, maxDose, num=maxDose+1)
0214 sum=np.array(sum)
0215 outputfile1.close()
0216
0217 print("\nCalculation concluded. Check the final image and close it to end the program.!")
0218
0219 plt.yscale('log')
0220 plt.xlim([0, 8])
0221 plt.ylim([1e-2, 1])
0222 plt.plot(c,np.exp(-sum[:]))
0223 plt.show()
0224
0225 print("\nThank you for using molecularDNAsurvival.!\n")
0226