| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Authors: O. Belov and M. Batmunkh
27 // January 2017
28 // last edit: L.T. Anh (2023)
29 /// \file BelovModel.cc
30 /// \brief Implementation of the BelovModel class
31
32 #include "BelovModel.hh"
33 #include "DamageClassifier.hh"
34 #include "ODESolver.hh"
35 #include <iostream>
36 #include <fstream>
37 #include <functional>
38 #include <limits>
39 #include <cmath>
40 #include <sstream>
41 #include <string>
42 #include <stdlib.h>
43 #include <time.h>
44 #include <ctime>
45 #include <vector>
46
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48
49 BelovModel::BelovModel():
50 fDz(0.),
51 falpha(0.),
52 fNirrep(0.),
53 fTime(0.)
54 {
55 for(int i=0;i<5;i++){
56 frepairsim[i].clear();
57 }
58 fdnarepair.clear();
59 fBpForDSB = 10;
60 }
61
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63
64 void BelovModel::Initialize()
65 {
66 for(int i=0;i<5;i++){
67 frepairsim[i].clear();
68 }
69 fdnarepair.clear();
70 }
71
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73
74 bool BelovModel::CalculateRepair(double Dz)
75 {
76 // Recalling model parameters
77 fDz = Dz;
78
79 std::cout << "Belov Model, CalculateRepair with:" << std::endl;
80 std::cout << " - Dz = " << fDz << " Gy" << std::endl;
81 std::cout << " - alpha = " << falpha << " Gy-1" << std::endl;
82 std::cout << " - Nirrep = " << fNirrep << std::endl;
83
84 if(falpha==0)
85 {
86 std::cout << " falpha=0:\n"
87 <<"- please use ComputeAndSetDamageInput function before calculating repair"
88 << std::endl
89 <<"- if above checked, then the reason might be: nDSBYiels is zero !!!"
90 << std::endl;
91 return false;
92 }
93
94 if(fNirrep==0)
95 {
96 std::cout << " fNirrep=0:\n"
97 <<"- please use ComputeAndSetDamageInput function before calculating repair"
98 << std::endl
99 <<"- if above checked, then the reason might be: nComplexDSBYiels is zero !!!"
100 << std::endl;
101 return false;
102 }
103
104 // INITIAL CONDITIONS
105
106 int NbEquat = 29; // Total number of model equations
107 std::vector<double> Y(NbEquat,0);
108
109 //---- Initial conditions for NHEJ -----
110
111 Y[0] = falpha;
112 Y[1] = Y[2] = Y[3] = Y[4] = Y[5] = Y[6] = Y[7] = Y[8] = Y[9] = 0.;
113
114 //---- Initial conditions for HR -------
115 Y[10] = Y[11] = Y[12] = Y[13] = Y[14] = Y[15] = Y[16] = Y[17]
116 = Y[18] = Y[19] = 0.;
117
118 //---- Initial conditions for SSA -----
119 Y[20] = Y[21] = Y[22] = Y[23] = Y[24] = 0.;
120
121 //---- Initial conditions for Alt-NHEJ (MMEJ) -----
122 Y[25] = Y[26] = Y[27] = Y[28] = 0.;
123
124 // Integration parameters
125 double t0 = 0.0; // Starting time point (dimensionless)
126 double t1 = 45.3; // Final time point (dimensionless)
127 double dt = 2.e-6; // Intergration time step (dimensionless)
128
129 double K8 = 0.552; // [h-1], scaling variable
130
131 std::function<std::vector<double>(double,std::vector<double>)>
132 func = [this] (double t,std::vector<double> y) -> std::vector<double> {
133 return Belov_odes_system(t,y);
134 };
135
136 std::vector<std::vector<double>> Y_vec;
137 std::vector<double> times;
138 double epsilon = 0.1;
139 ODESolver odeSolver;
140 odeSolver.SetNstepsForObserver(16000);
141 odeSolver.Embedded_RungeKutta_Fehlberg(func,Y,t0,t1,dt,epsilon,×,&Y_vec);
142 //odeSolver.RungeKutta4(func,Y,t0,t1,dt,×,&Y_vec);
143 size_t steps = Y_vec.size();
144
145 // Output options for different repair stages
146 size_t Nfoci=5;
147 std::string FociName;
148 double maxY= -1e-9;
149
150 for (size_t ifoci=0;ifoci<Nfoci;ifoci++)
151 {
152 maxY = -1e-9;
153
154 if(ifoci==0)FociName = std::string("Ku" );
155 if(ifoci==1)FociName = std::string("DNAPKcs");
156 if(ifoci==2)FociName = std::string("RPA" );
157 if(ifoci==3)FociName = std::string("Rad51" );
158 if(ifoci==4)FociName = std::string("gH2AX" );
159 for( size_t istep=0; istep<steps; istep+=1 )
160 {
161 double val = 0;
162 if(FociName=="Ku" ) val = Y_vec[istep][1];
163 if(FociName=="DNAPKcs") {
164 val = Y_vec[istep][3] +Y_vec[istep][4] +Y_vec[istep][5]+Y_vec[istep][6]+Y_vec[istep][7];
165 }
166 if(FociName=="RPA" ) val = Y_vec[istep][14]+Y_vec[istep][15]+Y_vec[istep][20];
167 if(FociName=="Rad51" ) val = Y_vec[istep][15]+Y_vec[istep][16]+Y_vec[istep][17];
168 if(FociName=="gH2AX" ) val = Y_vec[istep][9];
169 if (maxY < val ) maxY = val;
170 double time = times[istep]*K8;
171 frepairsim[ifoci].push_back(std::make_pair(time,val));
172 }
173 fdnarepair.insert(make_pair(FociName,frepairsim[ifoci]));
174 }
175 Y_vec.clear();Y_vec.shrink_to_fit();
176 return true;
177 }
178
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180
181
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
183
184 void BelovModel::ComputeAndSetDamageInput(std::vector<Damage> vecDamage)
185 {
186
187 DamageClassifier damClass;
188 auto classifiedDamage = damClass.MakeCluster(vecDamage,fBpForDSB,false);
189
190 ComplexDSBYield = damClass.GetNumComplexDSB(classifiedDamage);
191 DSBYield = damClass.GetNumDSB(classifiedDamage);
192 falpha = DSBYield/fDose;
193
194 fNirrep = ComplexDSBYield/DSBYield;
195 }
196
197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
198
199 std::vector<double> BelovModel::Belov_odes_system(double t,std::vector<double> Y)
200 {
201 // DSBRepairPathways
202 std::vector<double> YP;
203 // Concentrations of repair enzymes set to be constant
204 double X1; // [Ku]
205 double X2; // [DNAPKcsArt]
206 double X3; // [LigIV/XRCC4/XLF]
207 double X4; // [PNKP]
208 double X5; // [Pol]
209 double X6; // [H2AX]
210 double X7; // [MRN/CtIP/ExoI/Dna2]
211 double X8; // [ATM]
212 double X9; // [RPA]
213 double X10; // [Rad51/Rad51par/BRCA2]
214 double X11; // [DNAinc]
215 double X12; // [Rad52]
216 double X13; // [ERCC1/XPF]
217 double X14; // [LigIII]
218 double X15; // [PARP1]
219 double X16; // [Pol]
220 double X17; // [LigI]
221
222 X1 = X2 = X3 = X4 = X5 = X6 = X7 = X8 =
223 X9 = X10 = X11 = X12 = X13 = X14 =
224 X15 = X16 = X17 = 400000.;
225
226 fTime = t; // Recalling t
227
228 // DIMENSIONAL REACTION RATES
229 //
230 //------------NHEJ--------------
231 double K1 = 11.052; // M-1*h-1
232 double Kmin1 = 6.59999*1e-04; // h-1
233 double K2 = 18.8305*(1.08517-std::exp(-21.418/std::pow(fDz,1.822))); // M-1*h-1
234 double Kmin2 = 5.26*1e-01; //h-1
235 double K3 = 1.86; // h-1
236 double K4 = 1.38*1e+06; // M-1*h-1
237 double Kmin4 = 3.86*1e-04; // h-1
238 double K5 = 15.24; // M-1*h-1
239 double Kmin5 = 8.28; // h-1
240 double K6 = 18.06; // M-1*h-1
241 double Kmin6 = 1.33; // h-1
242 double K7 = 2.73*1e+05; // M-1*h-1
243 double Kmin7 = 3.2; // h-1
244 double K8 = 5.52*1e-01; // h-1
245 double K9 = 1.66*1e-01; // h-1
246 double K10 = (1.93*1e-07)/fNirrep; // M
247 double K11 = 7.50*1e-02; // h-1
248 double K12 = 11.1; // h-1
249 //
250 //------------HR--------------
251 double P1 = 1.75*1e+03; // M-1*h-1
252 double Pmin1 = 1.33*1e-04; // h-1
253 double P2 = 0.39192; // h-1
254 double Pmin2 = 2.7605512*1e+02; // h-1
255 double P3 = 1.37*1e+04; // M-1*h-1
256 double Pmin3 = 2.34; // h-1
257 double P4 = 3.588*1e-02; // h-1
258 double P5 = 1.20*1e+05; // M-1*h-1
259 double Pmin5 = 8.82*1e-05; // h-1
260 double P6 = 1.54368*1e+06; // M-1*h-1
261 double Pmin6 = 1.55*1e-03; // h-1
262 double P7 = 1.4904; // h-1
263 double P8 = 1.20*1e+04; // M-1*h-1
264 double Pmin8 = 2.49*1e-04; // h-1
265 double P9 = 1.104; //h-1
266 double P10 = 7.20*1e-03; // h-1
267 double P11 = 6.06*1e-04; // h-1
268 double P12 = 2.76*1e-01; // h-1
269 //
270 //------------SSA--------------
271 double Q1 = 1.9941*1e+05; // M-1*h-1
272 double Qmin1 = 1.71*1e-04; // h-1
273 double Q2 = 4.8052*1e+04; // M-1*h-1
274 double Q3 = 6*1e+03; // M-1*h-1
275 double Qmin3 = 6.06*1e-04; // h-1
276 double Q4 = 1.62*1e-03; // h-1
277 double Q5 = 8.40*1e+04; // M-1*h-1
278 double Qmin5 = 4.75*1e-04; // h-1
279 double Q6 = 11.58; // h-1
280 //
281 //-------alt-NHEJ (MMEJ)--------
282 double R1 = 2.39*1e+03; // M-1*h-1
283 double Rmin1 = 12.63; // h-1
284 double R2 = 4.07*1e+04; // M-1*h-1
285 double R3 = 9.82; // h-1
286 double R4 = 1.47*1e+05; // M-1*h-1
287 double Rmin4 = 2.72; // h-1
288 double R5 = 1.65*1e-01; //h-1
289 //
290 // Scalling rate XX1
291 double XX1 = 9.19*1e-07; // M
292 //
293 // DIMENSIONLESS REACTION RATES
294 //
295 //------------NHEJ--------------
296 double k1 = K1*XX1/K8;
297 double kmin1 = Kmin1/K8;
298 double k2 = K2*XX1/K8;
299 double kmin2 = Kmin2/K8;
300 double k3 = K3/K8;
301 double k4 = K4*XX1/K8;
302 double kmin4 = Kmin4/K8;
303 double k5 = K5*XX1/K8;
304 double kmin5 = Kmin5/K8;
305 double k6 = K6*XX1/K8;
306 double kmin6 = Kmin6/K8;
307 double k7 = K7*XX1/K8;
308 double kmin7 = Kmin7/K8;
309 double k8 = K8/K8;
310 double k9 = K9/K8;
311 double k10 = K10/XX1;
312 double k11 = K11/K8;
313 double k12 = K12/K8;
314 //
315 //------------HR--------------
316 double p1 = P1*XX1/K8;
317 double pmin1 = Pmin1/K8;
318 double p2 = P2/K8;
319 double pmin2 = Pmin2/K8;
320 double p3 = P3*XX1/K8;
321 double pmin3 = Pmin3/K8;
322 double p4 = P4/K8;
323 double p5 = P5*XX1/K8;
324 double pmin5 = Pmin5/K8;
325 double p6 = P6*XX1/K8;
326 double pmin6 = Pmin6/K8;
327 double p7 = P7/K8;
328 double p8 = P8*XX1/K8;
329 double pmin8 = Pmin8/K8;
330 double p9 = P9/K8;
331 double p10 = P10/K8;
332 double p11 = P11/K8;
333 double p12 = P12/K8;
334 //
335 //------------SSA--------------
336 double q1= Q1*XX1/K8;
337 double qmin1 = Qmin1/K8;
338 double q2 = Q2*XX1/K8;
339 double q3 = Q3*XX1/K8;
340 double qmin3 = Qmin3/K8;
341 double q4 = Q4/K8;
342 double q5 = Q5*XX1/K8;
343 double qmin5 = Qmin5/K8;
344 double q6 = Q6/K8;
345 //
346 //-------alt-NHEJ (MMEJ)--------
347 double r1 = R1*XX1/K8;
348 double rmin1 = Rmin1/K8;
349 double r2 = R2*XX1/K8;
350 double r3 = R3/K8;
351 double r4 = R4*XX1/K8;
352 double rmin4 = Rmin4/K8;
353 double r5 = R5/K8;
354 //------------------------------------
355
356 // SYSTEM OF DIFFERENTIAL EQUATIONS
357
358 // ----- NHEJ ----------
359 YP.push_back( fNirrep - k1*Y[0]*X1 + kmin1*Y[1] - p1*Y[0]*X1 + pmin1*Y[10]); // [DSB]
360
361 YP.push_back( k1*Y[0]*X1 - kmin1*Y[1] - k2*Y[1]*X2 + kmin2*Y[2]); // [DBS * Ku]
362
363 YP.push_back( k2*Y[1]*X2 - k3*Y[2] - kmin2*Y[2]); // [DSB * DNA-PK/Art]
364
365 YP.push_back( k3*Y[2] - k4*(Y[3]*Y[3]) + kmin4*Y[4]); // [DSB * DNA-PK/ArtP]
366
367 YP.push_back( k4*(Y[3]*Y[3]) - kmin4*Y[4] - k5*Y[4]*X3 + kmin5*Y[5]); // [Bridge]
368
369 YP.push_back( kmin6*Y[6] + k5*Y[4]*X3 - kmin5*Y[5] - k6*Y[5]*X4);
370 // [Bridge * LigIV/XRCC4/XLF]
371
372 YP.push_back( -kmin6*Y[6] - k7*Y[6]*X5 + kmin7*Y[7] + k6*Y[5]*X4);
373 // [Bridge * LigIV/XRCC4/XLF * PNKP]
374
375 YP.push_back( k7*Y[6]*X5 - k8*Y[7] - kmin7*Y[7]);
376 // [Bridge * LigIV/XRCC4/XLF * PNKP * Pol]
377
378 YP.push_back( r5*Y[28] + k8*Y[7] + p12*Y[18] + p11*Y[19] + q6*Y[24]); // [dsDNA]
379
380 YP.push_back( (k9*(Y[3] + Y[4] + Y[5] + Y[6] + Y[7])*X6)/(k10 + Y[3] + Y[4] + Y[5]
381 + Y[6] + Y[7]) - k11*Y[8] - k12*Y[9]); // [gH2AX foci]
382
383 // ----- HR ----------
384 YP.push_back( p1*Y[0]*X7 - pmin1*Y[10] - p3*Y[10]*Y[11] + pmin3*Y[12]);
385 // [MRN/CtIP/ExoI/Dna2]
386
387 YP.push_back( p2*X8 - pmin2*Y[11] - p3*Y[10]*Y[11] + p4*Y[12] + pmin3*Y[12]);
388 // [ATMP]
389
390 YP.push_back( p3*Y[10]*Y[11] - p4*Y[12] - pmin3*Y[12]);
391 // [DSB * MRN/CtIP/ExoI/Dna2 * ATMP]
392
393 YP.push_back( rmin1*Y[25] + p4*Y[12] - r1*X15*Y[13] - p5*Y[13]*X9 + pmin5*Y[14]);
394 // [ssDNA]
395
396 YP.push_back( pmin6*Y[15] + p5*Y[13]*X9 - pmin5*Y[14] - p6*Y[14]*X10 -
397 q1*Y[14]*X12 + qmin1*Y[20]); // [ssDNA * RPA]
398
399 YP.push_back( -p7*Y[15] - pmin6*Y[15] + p6*Y[14]*X10);
400 // [ssDNA * RPA * Rad51/Rad51par/BRCA2]
401
402 YP.push_back( p7*Y[15] - p8*Y[16]*X11 + pmin8*Y[17]); // [Rad51 filament]
403
404 YP.push_back( p8*Y[16]*X11 - p9*Y[17] - pmin8*Y[17]); // [Rad51 filament * DNAinc]
405
406 YP.push_back( p9*Y[17] - p10*Y[18] - p12*Y[18]); // [D-loop]
407
408 YP.push_back( p10*Y[18] - p11*Y[19]); // [dHJ]
409
410 // ----- SSA ----------
411 YP.push_back( q1*Y[14]*X12 - qmin1*Y[20] - q2*(Y[20]*Y[20]));
412 // [ssDNA * RPA * Rad52]
413
414 YP.push_back( q2*(Y[20]*Y[20]) - q3*Y[21]*X13 + qmin3*Y[22]); // [Flap]
415
416 YP.push_back( q3*Y[21]*X13 - q4*Y[22] - qmin3*Y[22]); // [Flap * ERCC1/XPF]
417
418 YP.push_back( q4*Y[22] - q5*Y[23]*X14 + qmin5*Y[24]); // [dsDNAnicks]
419
420 YP.push_back( q5*Y[23]*X14 - q6*Y[24] - qmin5*Y[24]); // [dsDNAnicks * LigIII]
421
422 // ----- MMEJ ----------
423 YP.push_back( -rmin1*Y[25] - r2*Y[25]*X16 + r1*X15*Y[13]); // [ssDNA * PARP1]
424
425 YP.push_back( r2*Y[25]*X16 - r3*Y[26]); // [ssDNA * Pol]
426
427 YP.push_back( r3*Y[26] - r4*Y[27]*X17 + rmin4*Y[28]); // [MicroHomol]
428
429 YP.push_back( r4*Y[27]*X17 - r5*Y[28] - rmin4*Y[28]); // [MicroHomol * LigI]
430
431 //---------------------------------------------------
432 return YP;
433 }
434
435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
436
437 std::vector<std::pair<double,double>> BelovModel::GetDNARepair(std::string NameFoci)
438 {
439 decltype(fdnarepair)::iterator it = fdnarepair.find(NameFoci);
440 if (it != fdnarepair.end()) {
441 return it->second;
442 }
443 else{
444 std::cerr<<"There is no Foci with name: "<<NameFoci<<" !!!"<<std::endl;
445 exit(0); //exception needed
446 }
447 }
448
449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
450
451 void BelovModel::WriteOutput(std::string pFileName)
452 {
453 std::fstream file;
454 file.open(pFileName.c_str(), std::ios_base::out);
455 //Header part
456 file <<"#===================================== BELOV MODEL ========================================#\n";
457 file << " Belov Model, CalculateRepair with:\n";
458 file << "#DSB = " << DSBYield << " (SB) " << "#Complex DSB= " << ComplexDSBYield << " (SB)\n";
459 file << "#Dz = " << fDz << " Gy\n";
460 file << "#Nirrep = " << fNirrep << "\n";
461 file <<"#===========================================================================================#\n";
462 file << "Time\t";
463 for(auto it=fdnarepair.begin();it!=fdnarepair.end();it++)
464 {
465 file << it->first << "\t";
466 }
467 file << "\n";
468 //End header part
469 int nVal = fdnarepair.begin()->second.size();
470
471 for(int i=0;i<nVal;i++)
472 {
473 file << fdnarepair["DNAPKcs"][i].first << "\t";
474
475 for(auto it=fdnarepair.begin();it!=fdnarepair.end();it++)
476 {
477 auto data = it->second;
478 file << data[i].second << "\t";
479 }
480 file << "\n";
481
482 }
483
484 file.close();
485 }
486
487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
488
489 void BelovModel::SetDSBandComDSBandDose(double dsb,double cdsb, double d)
490 {
491 SetDose(d);
492 ComplexDSBYield = cdsb;
493 DSBYield = dsb;
494 falpha = DSBYield/fDose;
495 fNirrep = ComplexDSBYield/DSBYield;
496 }
497
498 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
499