repairmodels/src/BelovModel.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Authors: O. Belov and M. Batmunkh
0027 // January 2017
0028 // last edit: L.T. Anh (2023)
0029 /// \file BelovModel.cc
0030 /// \brief Implementation of the BelovModel class
0031
0032 #include "BelovModel.hh"
0033 #include "DamageClassifier.hh"
0034 #include "ODESolver.hh"
0035 #include <iostream>
0036 #include <fstream>
0037 #include <functional>
0038 #include <limits>
0039 #include <cmath>
0040 #include <sstream>
0041 #include <string>
0042 #include <stdlib.h>
0043 #include <time.h>
0044 #include <ctime>
0045 #include <vector>
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048
0049 BelovModel::BelovModel():
0050 fDz(0.),
0051 falpha(0.),
0052 fNirrep(0.),
0053 fTime(0.)
0054 {
0055   for(int i=0;i<5;i++){
0056     frepairsim[i].clear();
0057   }
0058   fdnarepair.clear();
0059   fBpForDSB = 10;
0060 }
0061
0062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0063
0064 void BelovModel::Initialize()
0065 {
0066   for(int i=0;i<5;i++){
0067     frepairsim[i].clear();
0068   }
0069   fdnarepair.clear();
0070 }
0071
0072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0073
0074 bool BelovModel::CalculateRepair(double Dz)
0075 {
0076     // Recalling model parameters
0077     fDz = Dz;
0078
0079     std::cout << "Belov Model, CalculateRepair with:" << std::endl;
0080     std::cout << "   - Dz = " << fDz << " Gy" << std::endl;
0081     std::cout << "   - alpha = " << falpha << " Gy-1" << std::endl;
0082     std::cout << "   - Nirrep = " << fNirrep << std::endl;
0083
0084     if(falpha==0)
0085     {
0086     std::cout << " falpha=0:\n"
0087                 <<"- please use ComputeAndSetDamageInput function before calculating repair"
0088                 << std::endl
0089                 <<"- if above checked, then the reason might be: nDSBYiels is zero !!!"
0090                 << std::endl;
0091     return false;
0092     }
0093
0094     if(fNirrep==0)
0095     {
0096     std::cout << " fNirrep=0:\n"
0097                 <<"- please use ComputeAndSetDamageInput function before calculating repair"
0098                 << std::endl
0099                 <<"- if above checked, then the reason might be: nComplexDSBYiels is zero !!!"
0100                 << std::endl;
0101     return false;
0102     }
0103
0104     // INITIAL CONDITIONS
0105
0106     int NbEquat = 29;   // Total number of model equations
0107     std::vector<double> Y(NbEquat,0);
0108
0109     //---- Initial conditions for NHEJ -----
0110
0111     Y[0] = falpha;
0112     Y[1] = Y[2] = Y[3] = Y[4] = Y[5] = Y[6] = Y[7] = Y[8] = Y[9] = 0.;
0113
0114     //---- Initial conditions for HR -------
0115     Y[10] = Y[11] = Y[12] = Y[13] = Y[14] = Y[15] = Y[16] = Y[17]
0116     = Y[18] = Y[19] = 0.;
0117
0118     //---- Initial conditions for SSA -----
0119     Y[20] = Y[21] = Y[22] = Y[23] = Y[24] = 0.;
0120
0121     //---- Initial conditions for Alt-NHEJ (MMEJ) -----
0122     Y[25] = Y[26] = Y[27] = Y[28] = 0.;
0123
0124     // Integration parameters
0125     double t0 = 0.0;   // Starting time point (dimensionless)
0126     double t1 = 45.3;  // Final time point (dimensionless)
0127     double dt = 2.e-6; // Intergration time step (dimensionless)
0128
0129     double K8 = 0.552; // [h-1], scaling variable
0130
0131     std::function<std::vector<double>(double,std::vector<double>)>
0132     func = [this] (double t,std::vector<double> y) -> std::vector<double> {
0133         return Belov_odes_system(t,y);
0134     };
0135
0136     std::vector<std::vector<double>> Y_vec;
0137     std::vector<double> times;
0138     double epsilon = 0.1;
0139     ODESolver odeSolver;
0140     odeSolver.SetNstepsForObserver(16000);
0141     odeSolver.Embedded_RungeKutta_Fehlberg(func,Y,t0,t1,dt,epsilon,&times,&Y_vec);
0142     //odeSolver.RungeKutta4(func,Y,t0,t1,dt,&times,&Y_vec);
0143     size_t steps = Y_vec.size();
0144
0145     // Output options for different repair stages
0146     size_t Nfoci=5;
0147     std::string FociName;
0148     double maxY= -1e-9;
0149
0150     for (size_t ifoci=0;ifoci<Nfoci;ifoci++)
0151     {
0152     maxY = -1e-9;
0153
0154     if(ifoci==0)FociName = std::string("Ku"     );
0155     if(ifoci==1)FociName = std::string("DNAPKcs");
0156     if(ifoci==2)FociName = std::string("RPA"    );
0157     if(ifoci==3)FociName = std::string("Rad51"  );
0158     if(ifoci==4)FociName = std::string("gH2AX"  );
0159     for( size_t istep=0; istep<steps; istep+=1 )
0160     {
0161         double val = 0;
0162         if(FociName=="Ku"     ) val = Y_vec[istep][1];
0163         if(FociName=="DNAPKcs") {
0164             val = Y_vec[istep][3] +Y_vec[istep][4] +Y_vec[istep][5]+Y_vec[istep][6]+Y_vec[istep][7];
0165         }
0166         if(FociName=="RPA"    ) val = Y_vec[istep][14]+Y_vec[istep][15]+Y_vec[istep][20];
0167         if(FociName=="Rad51"  ) val = Y_vec[istep][15]+Y_vec[istep][16]+Y_vec[istep][17];
0168         if(FociName=="gH2AX"  ) val = Y_vec[istep][9];
0169         if (maxY < val ) maxY = val;
0170         double time = times[istep]*K8;
0171         frepairsim[ifoci].push_back(std::make_pair(time,val));
0172     }
0173     fdnarepair.insert(make_pair(FociName,frepairsim[ifoci]));
0174     }
0175     Y_vec.clear();Y_vec.shrink_to_fit();
0176     return true;
0177 }
0178
0179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0180
0181
0182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0183
0184 void BelovModel::ComputeAndSetDamageInput(std::vector<Damage> vecDamage)
0185 {
0186
0187   DamageClassifier damClass;
0188   auto classifiedDamage = damClass.MakeCluster(vecDamage,fBpForDSB,false);
0189
0190   ComplexDSBYield = damClass.GetNumComplexDSB(classifiedDamage);
0191   DSBYield = damClass.GetNumDSB(classifiedDamage);
0192   falpha = DSBYield/fDose;
0193
0194   fNirrep = ComplexDSBYield/DSBYield;
0195 }
0196
0197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0198
0199 std::vector<double> BelovModel::Belov_odes_system(double t,std::vector<double> Y)
0200 {
0201     // DSBRepairPathways
0202     std::vector<double> YP;
0203     // Concentrations of repair enzymes set to be constant
0204     double X1; // [Ku]
0205     double X2; // [DNAPKcsArt]
0206     double X3; // [LigIV/XRCC4/XLF]
0207     double X4; // [PNKP]
0208     double X5; // [Pol]
0209     double X6; // [H2AX]
0210     double X7; // [MRN/CtIP/ExoI/Dna2]
0211     double X8; // [ATM]
0212     double X9; // [RPA]
0213     double X10; // [Rad51/Rad51par/BRCA2]
0214     double X11; // [DNAinc]
0215     double X12; // [Rad52]
0216     double X13; // [ERCC1/XPF]
0217     double X14; // [LigIII]
0218     double X15; // [PARP1]
0219     double X16; // [Pol]
0220     double X17; // [LigI]
0221
0222     X1 = X2 = X3 = X4 = X5 = X6 = X7 = X8 =
0223     X9 = X10 = X11 = X12 = X13 = X14 =
0224     X15 = X16 = X17 = 400000.;
0225
0226     fTime = t;  // Recalling t
0227
0228     //  DIMENSIONAL REACTION RATES
0229     //
0230     //------------NHEJ--------------
0231     double K1 = 11.052;             // M-1*h-1
0232     double Kmin1 = 6.59999*1e-04;   // h-1
0233     double K2 = 18.8305*(1.08517-std::exp(-21.418/std::pow(fDz,1.822))); // M-1*h-1
0234     double Kmin2 = 5.26*1e-01;      //h-1
0235     double K3 = 1.86;               // h-1
0236     double K4 = 1.38*1e+06;          // M-1*h-1
0237     double Kmin4 = 3.86*1e-04;      // h-1
0238     double K5 = 15.24;              // M-1*h-1
0239     double Kmin5 = 8.28;            // h-1
0240     double K6 = 18.06;              // M-1*h-1
0241     double Kmin6 = 1.33;            // h-1
0242     double K7 = 2.73*1e+05;          // M-1*h-1
0243     double Kmin7 = 3.2;             // h-1
0244     double K8 = 5.52*1e-01;         // h-1
0245     double K9 = 1.66*1e-01;         // h-1
0246     double K10 = (1.93*1e-07)/fNirrep; // M
0247     double K11 = 7.50*1e-02;        // h-1
0248     double K12 = 11.1;              // h-1
0249     //
0250     //------------HR--------------
0251     double P1 = 1.75*1e+03;          // M-1*h-1
0252     double Pmin1 = 1.33*1e-04;      // h-1
0253     double P2 = 0.39192;            // h-1
0254     double Pmin2 = 2.7605512*1e+02;  // h-1
0255     double P3 = 1.37*1e+04;          // M-1*h-1
0256     double Pmin3 = 2.34;            // h-1
0257     double P4 = 3.588*1e-02;       // h-1
0258     double P5 = 1.20*1e+05;          // M-1*h-1
0259     double Pmin5 = 8.82*1e-05;      // h-1
0260     double P6 = 1.54368*1e+06;       // M-1*h-1
0261     double Pmin6 = 1.55*1e-03;      // h-1
0262     double P7 = 1.4904;              // h-1
0263     double P8 = 1.20*1e+04;          // M-1*h-1
0264     double Pmin8 = 2.49*1e-04;      // h-1
0265     double P9 = 1.104;            //h-1
0266     double P10 = 7.20*1e-03;        // h-1
0267     double P11 = 6.06*1e-04;        // h-1
0268     double P12 = 2.76*1e-01;        // h-1
0269     //
0270     //------------SSA--------------
0271     double Q1 = 1.9941*1e+05;       // M-1*h-1
0272     double Qmin1 = 1.71*1e-04;      // h-1
0273     double Q2 = 4.8052*1e+04;       // M-1*h-1
0274     double Q3 = 6*1e+03;             // M-1*h-1
0275     double Qmin3 = 6.06*1e-04;      // h-1
0276     double Q4 = 1.62*1e-03;         // h-1
0277     double Q5 = 8.40*1e+04;          // M-1*h-1
0278     double Qmin5 = 4.75*1e-04;      // h-1
0279     double Q6 = 11.58;              // h-1
0280     //
0281     //-------alt-NHEJ (MMEJ)--------
0282     double R1 = 2.39*1e+03;          // M-1*h-1
0283     double Rmin1 = 12.63;           // h-1
0284     double R2 = 4.07*1e+04;          // M-1*h-1
0285     double R3 = 9.82;               // h-1
0286     double R4 = 1.47*1e+05;          // M-1*h-1
0287     double Rmin4 = 2.72;            // h-1
0288     double R5 = 1.65*1e-01;         //h-1
0289     //
0290     // Scalling rate XX1
0291     double XX1 = 9.19*1e-07; // M
0292     //
0293     // DIMENSIONLESS REACTION RATES
0294     //
0295     //------------NHEJ--------------
0296     double k1 = K1*XX1/K8;
0297     double kmin1 = Kmin1/K8;
0298     double k2 = K2*XX1/K8;
0299     double kmin2 = Kmin2/K8;
0300     double k3 = K3/K8;
0301     double k4 = K4*XX1/K8;
0302     double kmin4 = Kmin4/K8;
0303     double k5 = K5*XX1/K8;
0304     double kmin5 = Kmin5/K8;
0305     double k6 = K6*XX1/K8;
0306     double kmin6 = Kmin6/K8;
0307     double k7 = K7*XX1/K8;
0308     double kmin7 = Kmin7/K8;
0309     double k8 = K8/K8;
0310     double k9 = K9/K8;
0311     double k10 = K10/XX1;
0312     double k11 = K11/K8;
0313     double k12 = K12/K8;
0314     //
0315     //------------HR--------------
0316     double p1 = P1*XX1/K8;
0317     double pmin1 = Pmin1/K8;
0318     double p2 = P2/K8;
0319     double pmin2 = Pmin2/K8;
0320     double p3 = P3*XX1/K8;
0321     double pmin3 = Pmin3/K8;
0322     double p4 = P4/K8;
0323     double p5 = P5*XX1/K8;
0324     double pmin5 = Pmin5/K8;
0325     double p6 = P6*XX1/K8;
0326     double pmin6 = Pmin6/K8;
0327     double p7 = P7/K8;
0328     double p8 = P8*XX1/K8;
0329     double pmin8 = Pmin8/K8;
0330     double p9 = P9/K8;
0331     double p10 = P10/K8;
0332     double p11 = P11/K8;
0333     double p12 = P12/K8;
0334     //
0335     //------------SSA--------------
0336     double q1= Q1*XX1/K8;
0337     double qmin1 = Qmin1/K8;
0338     double q2 = Q2*XX1/K8;
0339     double q3 = Q3*XX1/K8;
0340     double qmin3 = Qmin3/K8;
0341     double q4 = Q4/K8;
0342     double q5 = Q5*XX1/K8;
0343     double qmin5 = Qmin5/K8;
0344     double q6 = Q6/K8;
0345     //
0346     //-------alt-NHEJ (MMEJ)--------
0347     double r1 = R1*XX1/K8;
0348     double rmin1 = Rmin1/K8;
0349     double r2 = R2*XX1/K8;
0350     double r3 = R3/K8;
0351     double r4 = R4*XX1/K8;
0352     double rmin4 = Rmin4/K8;
0353     double r5 = R5/K8;
0354     //------------------------------------
0355
0356     // SYSTEM OF DIFFERENTIAL EQUATIONS
0357
0358     // ----- NHEJ ----------
0359     YP.push_back( fNirrep - k1*Y[0]*X1 + kmin1*Y[1] - p1*Y[0]*X1 + pmin1*Y[10]); // [DSB]
0360
0361     YP.push_back( k1*Y[0]*X1 - kmin1*Y[1] - k2*Y[1]*X2 + kmin2*Y[2]); // [DBS * Ku]
0362
0363     YP.push_back( k2*Y[1]*X2 - k3*Y[2] - kmin2*Y[2]); // [DSB * DNA-PK/Art]
0364
0365     YP.push_back( k3*Y[2] - k4*(Y[3]*Y[3]) + kmin4*Y[4]); // [DSB * DNA-PK/ArtP]
0366
0367     YP.push_back( k4*(Y[3]*Y[3]) - kmin4*Y[4] - k5*Y[4]*X3 + kmin5*Y[5]); // [Bridge]
0368
0369     YP.push_back( kmin6*Y[6] + k5*Y[4]*X3 - kmin5*Y[5] - k6*Y[5]*X4);
0370     // [Bridge * LigIV/XRCC4/XLF]
0371
0372     YP.push_back( -kmin6*Y[6] - k7*Y[6]*X5 +  kmin7*Y[7] + k6*Y[5]*X4);
0373     // [Bridge * LigIV/XRCC4/XLF * PNKP]
0374
0375     YP.push_back( k7*Y[6]*X5 - k8*Y[7] - kmin7*Y[7]);
0376     // [Bridge * LigIV/XRCC4/XLF * PNKP * Pol]
0377
0378     YP.push_back( r5*Y[28] + k8*Y[7] + p12*Y[18] + p11*Y[19] + q6*Y[24]); // [dsDNA]
0379
0380     YP.push_back( (k9*(Y[3] + Y[4] + Y[5] + Y[6] + Y[7])*X6)/(k10 + Y[3] + Y[4] + Y[5]
0381             + Y[6] + Y[7]) - k11*Y[8] - k12*Y[9]); // [gH2AX foci]
0382
0383     // ----- HR ----------
0384     YP.push_back( p1*Y[0]*X7 - pmin1*Y[10] - p3*Y[10]*Y[11] + pmin3*Y[12]);
0385     // [MRN/CtIP/ExoI/Dna2]
0386
0387     YP.push_back( p2*X8 - pmin2*Y[11] - p3*Y[10]*Y[11] + p4*Y[12] + pmin3*Y[12]);
0388     // [ATMP]
0389
0390     YP.push_back( p3*Y[10]*Y[11] - p4*Y[12] - pmin3*Y[12]);
0391     // [DSB * MRN/CtIP/ExoI/Dna2 * ATMP]
0392
0393     YP.push_back( rmin1*Y[25] + p4*Y[12] - r1*X15*Y[13] - p5*Y[13]*X9 + pmin5*Y[14]);
0394     // [ssDNA]
0395
0396     YP.push_back( pmin6*Y[15] + p5*Y[13]*X9 - pmin5*Y[14] - p6*Y[14]*X10 -
0397             q1*Y[14]*X12 + qmin1*Y[20]); // [ssDNA * RPA]
0398
0399     YP.push_back( -p7*Y[15] - pmin6*Y[15] + p6*Y[14]*X10);
0400     // [ssDNA * RPA * Rad51/Rad51par/BRCA2]
0401
0402     YP.push_back( p7*Y[15] - p8*Y[16]*X11 + pmin8*Y[17]); // [Rad51 filament]
0403
0404     YP.push_back( p8*Y[16]*X11 - p9*Y[17] - pmin8*Y[17]); // [Rad51 filament * DNAinc]
0405
0406     YP.push_back( p9*Y[17] - p10*Y[18] - p12*Y[18]); // [D-loop]
0407
0408     YP.push_back( p10*Y[18] - p11*Y[19]); // [dHJ]
0409
0410     // ----- SSA ----------
0411     YP.push_back( q1*Y[14]*X12 - qmin1*Y[20] -  q2*(Y[20]*Y[20]));
0412     // [ssDNA * RPA * Rad52]
0413
0414     YP.push_back( q2*(Y[20]*Y[20]) - q3*Y[21]*X13 + qmin3*Y[22]); // [Flap]
0415
0416     YP.push_back( q3*Y[21]*X13 - q4*Y[22] - qmin3*Y[22]); // [Flap * ERCC1/XPF]
0417
0418     YP.push_back( q4*Y[22] - q5*Y[23]*X14 + qmin5*Y[24]); // [dsDNAnicks]
0419
0420     YP.push_back( q5*Y[23]*X14 - q6*Y[24] - qmin5*Y[24]); // [dsDNAnicks * LigIII]
0421
0422     // ----- MMEJ ----------
0423     YP.push_back( -rmin1*Y[25] - r2*Y[25]*X16 + r1*X15*Y[13]); // [ssDNA * PARP1]
0424
0425     YP.push_back( r2*Y[25]*X16 - r3*Y[26]); // [ssDNA * Pol]
0426
0427     YP.push_back( r3*Y[26] - r4*Y[27]*X17 + rmin4*Y[28]); // [MicroHomol]
0428
0429     YP.push_back( r4*Y[27]*X17 - r5*Y[28] - rmin4*Y[28]); // [MicroHomol * LigI]
0430
0431     //---------------------------------------------------
0432     return YP;
0433 }
0434
0435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0436
0437 std::vector<std::pair<double,double>> BelovModel::GetDNARepair(std::string NameFoci)
0438 {
0439     decltype(fdnarepair)::iterator it = fdnarepair.find(NameFoci);
0440   if (it != fdnarepair.end()) {
0441     return  it->second;
0442   }
0443   else{
0444     std::cerr<<"There is no Foci with name: "<<NameFoci<<" !!!"<<std::endl;
0445     exit(0);            //exception needed
0446   }
0447 }
0448
0449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0450
0451 void BelovModel::WriteOutput(std::string pFileName)
0452 {
0453   std::fstream file;
0454   file.open(pFileName.c_str(), std::ios_base::out);
0455   //Header part
0456   file <<"#===================================== BELOV  MODEL ========================================#\n";
0457   file << "                               Belov Model, CalculateRepair with:\n";
0458   file << "#DSB  = " <<  DSBYield   << " (SB)        " << "#Complex DSB= " << ComplexDSBYield << " (SB)\n";
0459   file << "#Dz     = " << fDz << " Gy\n";
0460   file << "#Nirrep = " << fNirrep << "\n";
0461   file <<"#===========================================================================================#\n";
0462   file << "Time\t";
0463   for(auto it=fdnarepair.begin();it!=fdnarepair.end();it++)
0464   {
0465     file << it->first << "\t";
0466   }
0467   file << "\n";
0468   //End header part
0469   int nVal = fdnarepair.begin()->second.size();
0470
0471   for(int i=0;i<nVal;i++)
0472   {
0473     file << fdnarepair["DNAPKcs"][i].first << "\t";
0474
0475     for(auto it=fdnarepair.begin();it!=fdnarepair.end();it++)
0476     {
0477       auto data = it->second;
0478       file << data[i].second << "\t";
0479     }
0480     file << "\n";
0481
0482   }
0483
0484   file.close();
0485 }
0486
0487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0488
0489 void BelovModel::SetDSBandComDSBandDose(double dsb,double cdsb, double d)
0490 {
0491     SetDose(d);
0492     ComplexDSBYield = cdsb;
0493     DSBYield = dsb;
0494     falpha = DSBYield/fDose;
0495     fNirrep = ComplexDSBYield/DSBYield;
0496 }
0497
0498 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......