| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 cl
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........oo
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........oo
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........oo
73
74 bool BelovModel::CalculateRepair(double Dz)
75 {
76 // Recalling model parameters
77 fDz = Dz;
78
79 std::cout << "Belov Model, CalculateRepair w
80 std::cout << " - Dz = " << fDz << " Gy" <<
81 std::cout << " - alpha = " << falpha << "
82 std::cout << " - Nirrep = " << fNirrep <<
83
84 if(falpha==0)
85 {
86 std::cout << " falpha=0:\n"
87 <<"- please use ComputeAndSetDamageInp
88 << std::endl
89 <<"- if above checked, then the reason
90 << std::endl;
91 return false;
92 }
93
94 if(fNirrep==0)
95 {
96 std::cout << " fNirrep=0:\n"
97 <<"- please use ComputeAndSetDamageInp
98 << std::endl
99 <<"- if above checked, then the reason
100 << std::endl;
101 return false;
102 }
103
104 // INITIAL CONDITIONS
105
106 int NbEquat = 29; // Total number of model
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[
113
114 //---- Initial conditions for HR -------
115 Y[10] = Y[11] = Y[12] = Y[13] = Y[14] = Y[15
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 (d
126 double t1 = 45.3; // Final time point (dime
127 double dt = 2.e-6; // Intergration time step
128
129 double K8 = 0.552; // [h-1], scaling variabl
130
131 std::function<std::vector<double>(double,std
132 func = [this] (double t,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,
142 //odeSolver.RungeKutta4(func,Y,t0,t1,dt,&tim
143 size_t steps = Y_vec.size();
144
145 // Output options for different repair stage
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]
163 if(FociName=="DNAPKcs") {
164 val = Y_vec[istep][3] +Y_vec[istep][4] +
165 }
166 if(FociName=="RPA" ) val = Y_vec[istep]
167 if(FociName=="Rad51" ) val = Y_vec[istep]
168 if(FociName=="gH2AX" ) val = Y_vec[istep]
169 if (maxY < val ) maxY = val;
170 double time = times[istep]*K8;
171 frepairsim[ifoci].push_back(std::make_pair
172 }
173 fdnarepair.insert(make_pair(FociName,frepair
174 }
175 Y_vec.clear();Y_vec.shrink_to_fit();
176 return true;
177 }
178
179 //....oooOO0OOooo........oooOO0OOooo........oo
180
181
182 //....oooOO0OOooo........oooOO0OOooo........oo
183
184 void BelovModel::ComputeAndSetDamageInput(std:
185 {
186
187 DamageClassifier damClass;
188 auto classifiedDamage = damClass.MakeCluster
189
190 ComplexDSBYield = damClass.GetNumComplexDSB(
191 DSBYield = damClass.GetNumDSB(classifiedDama
192 falpha = DSBYield/fDose;
193
194 fNirrep = ComplexDSBYield/DSBYield;
195 }
196
197 //....oooOO0OOooo........oooOO0OOooo........oo
198
199 std::vector<double> BelovModel::Belov_odes_sys
200 {
201 // DSBRepairPathways
202 std::vector<double> YP;
203 // Concentrations of repair enzymes set to b
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.41
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
360
361 YP.push_back( k1*Y[0]*X1 - kmin1*Y[1] - k2*Y
362
363 YP.push_back( k2*Y[1]*X2 - k3*Y[2] - kmin2*Y
364
365 YP.push_back( k3*Y[2] - k4*(Y[3]*Y[3]) + kmi
366
367 YP.push_back( k4*(Y[3]*Y[3]) - kmin4*Y[4] -
368
369 YP.push_back( kmin6*Y[6] + k5*Y[4]*X3 - kmin
370 // [Bridge * LigIV/XRCC4/XLF]
371
372 YP.push_back( -kmin6*Y[6] - k7*Y[6]*X5 + km
373 // [Bridge * LigIV/XRCC4/XLF * PNKP]
374
375 YP.push_back( k7*Y[6]*X5 - k8*Y[7] - kmin7*Y
376 // [Bridge * LigIV/XRCC4/XLF * PNKP * Pol]
377
378 YP.push_back( r5*Y[28] + k8*Y[7] + p12*Y[18]
379
380 YP.push_back( (k9*(Y[3] + Y[4] + Y[5] + Y[6]
381 + Y[6] + Y[7]) - k11*Y[8] - k12*Y[9]); /
382
383 // ----- HR ----------
384 YP.push_back( p1*Y[0]*X7 - pmin1*Y[10] - p3*
385 // [MRN/CtIP/ExoI/Dna2]
386
387 YP.push_back( p2*X8 - pmin2*Y[11] - p3*Y[10]
388 // [ATMP]
389
390 YP.push_back( p3*Y[10]*Y[11] - p4*Y[12] - pm
391 // [DSB * MRN/CtIP/ExoI/Dna2 * ATMP]
392
393 YP.push_back( rmin1*Y[25] + p4*Y[12] - r1*X1
394 // [ssDNA]
395
396 YP.push_back( pmin6*Y[15] + p5*Y[13]*X9 - pm
397 q1*Y[14]*X12 + qmin1*Y[20]); // [ssDNA *
398
399 YP.push_back( -p7*Y[15] - pmin6*Y[15] + p6*Y
400 // [ssDNA * RPA * Rad51/Rad51par/BRCA2]
401
402 YP.push_back( p7*Y[15] - p8*Y[16]*X11 + pmin
403
404 YP.push_back( p8*Y[16]*X11 - p9*Y[17] - pmin
405
406 YP.push_back( p9*Y[17] - p10*Y[18] - p12*Y[1
407
408 YP.push_back( p10*Y[18] - p11*Y[19]); // [dH
409
410 // ----- SSA ----------
411 YP.push_back( q1*Y[14]*X12 - qmin1*Y[20] -
412 // [ssDNA * RPA * Rad52]
413
414 YP.push_back( q2*(Y[20]*Y[20]) - q3*Y[21]*X1
415
416 YP.push_back( q3*Y[21]*X13 - q4*Y[22] - qmin
417
418 YP.push_back( q4*Y[22] - q5*Y[23]*X14 + qmin
419
420 YP.push_back( q5*Y[23]*X14 - q6*Y[24] - qmin
421
422 // ----- MMEJ ----------
423 YP.push_back( -rmin1*Y[25] - r2*Y[25]*X16 +
424
425 YP.push_back( r2*Y[25]*X16 - r3*Y[26]); // [
426
427 YP.push_back( r3*Y[26] - r4*Y[27]*X17 + rmin
428
429 YP.push_back( r4*Y[27]*X17 - r5*Y[28] - rmin
430
431 //------------------------------------------
432 return YP;
433 }
434
435 //....oooOO0OOooo........oooOO0OOooo........oo
436
437 std::vector<std::pair<double,double>> BelovMod
438 {
439 decltype(fdnarepair)::iterator it = fdnarepa
440 if (it != fdnarepair.end()) {
441 return it->second;
442 }
443 else{
444 std::cerr<<"There is no Foci with name: "<<N
445 exit(0); //exception needed
446 }
447 }
448
449 //....oooOO0OOooo........oooOO0OOooo........oo
450
451 void BelovModel::WriteOutput(std::string pFile
452 {
453 std::fstream file;
454 file.open(pFileName.c_str(), std::ios_base::
455 //Header part
456 file <<"#===================================
457 file << " Belo
458 file << "#DSB = " << DSBYield << " (SB)
459 file << "#Dz = " << fDz << " Gy\n";
460 file << "#Nirrep = " << fNirrep << "\n";
461 file <<"#===================================
462 file << "Time\t";
463 for(auto it=fdnarepair.begin();it!=fdnarepai
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 <<
474
475 for(auto it=fdnarepair.begin();it!=fdnarep
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........oo
488
489 void BelovModel::SetDSBandComDSBandDose(double
490 {
491 SetDose(d);
492 ComplexDSBYield = cdsb;
493 DSBYield = dsb;
494 falpha = DSBYield/fDose;
495 fNirrep = ComplexDSBYield/DSBYield;
496 }
497
498 //....oooOO0OOooo........oooOO0OOooo........oo
499