Skip to content

Extend IRT chemistry simulations for PBS buffer and additional salts #121

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
JoseRamosMendez merged 2 commits into from
Mar 18, 2026
Merged

Extend IRT chemistry simulations for PBS buffer and additional salts #121

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e2b2001
IRT Config: Implementation PBS and salt effects
Mar 17, 2026
f9d2400
IRT Config: PBS buffer components added and bugfixes
MarcBHahn Mar 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
108 changes: 103 additions & 5 deletions processes/TsIRTConfiguration.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -361,6 +361,22 @@ fKick(false), fAllTotallyDiffusionControlled(false)
if ( fPm->ParameterExists("Ch/"+chemistryList+"/ModelAcidPropertiesFromSubstance") ) {
fpHSolventConcentration = 0.0;
fpHValue = 7.1;

if (fPm->ParameterExists("Ch/"+chemistryList+"/MonovalentSalt")){
fMonovalentSalt = fPm->GetDoubleParameter("Ch/"+chemistryList+"/MonovalentSalt","molar concentration");
}
else{fMonovalentSalt=0.0;}

if (fPm->ParameterExists("Ch/"+chemistryList+"/DivalentSalt")){
fDivalentSalt = fPm->GetDoubleParameter("Ch/"+chemistryList+"/DivalentSalt","molar concentration");
}
else{fDivalentSalt=0.0;}

if (fPm->ParameterExists("Ch/"+chemistryList+"/TrivalentSalt")){
fTrivalentSalt = fPm->GetDoubleParameter("Ch/"+chemistryList+"/TrivalentSalt","molar concentration");
}
else{fTrivalentSalt=0.0;}


fpHSolvent = fPm->GetStringParameter("Ch/"+chemistryList+"/ModelAcidPropertiesFromSubstance");
G4StrUtil::to_lower(fpHSolvent);
Expand All @@ -380,6 +396,7 @@ fKick(false), fAllTotallyDiffusionControlled(false)
AdjustReactionRateForPH("PH");
}


G4cout << "======================================== performed ph Scaling ===========================================" << G4endl;
ResolveReactionParameters();
PrintReactionsInformation();
Expand Down Expand Up @@ -1023,6 +1040,37 @@ std::vector<G4double> TsIRTConfiguration::GetH2SO4ComponentsConcentrationPH(G4do
}


std::vector<G4double> TsIRTConfiguration::GetPhosphateBufferComponentsConcentrationPHandIonicStrength(G4double pH) {

G4double _Kw = 1E-14;
G4double _H_pos = pow(10,-pH);
G4double _OH_me = _Kw / _H_pos;

// only valid around pH 7.4 with 12mM Phosphate buffer component which is fixed for standard PBS
// pKa2 H_2PO_4^- <-> HPO_4^2- is 7.2 and the other two pKas are 2.15 and 12.35
// keep these fix, but include already for future extensions
G4double _H3PO_4 = 0.0;
G4double _PO_4 = 0.0;
G4double pKa1 = 2.1;
G4double pKa3 = 12.3;

// -> approximate by the Henderson Hasselbach equation pH=pKa + log_10(HPO_4^2-/H_2PO_4^-) (https://goldbook.iupac.org/terms/view/H02781)
// Reminder: adjust Ionic strength later on for standard PBS salts: 137 mM NaCl and 2.7mM KCl
G4double pKa2 = 7.2;
G4double buffer_component_concentration = 0.012; // M
G4double component_ratio = pow(10, pH-pKa2);
G4double _H_2PO_4 = buffer_component_concentration / (1.0 + component_ratio );
G4double _HPO_4 = component_ratio * _H_2PO_4;

G4double IonicStrength = 0.5 * (_H_pos + _H_2PO_4 + ( _HPO_4* 4) + ( _PO_4* 9) + _OH_me);
std::vector<G4double> Results = {IonicStrength, _H_pos , _H3PO_4, _H_2PO_4 , _HPO_4, _PO_4, _OH_me};
return Results;

}




std::vector<G4double> TsIRTConfiguration::GetH2SO4ComponentsConcentrationP(G4double Concentration) {
G4double Ka1 = pow(10,3);
G4double Ka2 = pow(10,-1.987);
Expand Down Expand Up @@ -1071,6 +1119,15 @@ G4double TsIRTConfiguration::GetIonicStrength(std::vector<G4double> Components)
return I;
}


G4double TsIRTConfiguration::GetIonicStrengthFromSalts(G4double fMonovalentSalt = 0.0, G4double fDivalentSalt = 0.0, G4double fTrivalentSalt = 0.0) {
//assuming all counterions with q=+-1
G4double counterions = fMonovalentSalt + 2.0 * fDivalentSalt + 3.0 * fTrivalentSalt;
G4double I = 0.5 * ( fMonovalentSalt + 4* fDivalentSalt + 9*fTrivalentSalt + counterions);
return I;
}


G4double TsIRTConfiguration::IonicRate(G4double IonicStrength, TsMolecularReaction Reaction) {
G4double a = 0.15;
G4int Charge1 = fMoleculesDefinition[Reaction.reactorA].charge;
Expand Down Expand Up @@ -1311,12 +1368,17 @@ void TsIRTConfiguration::AdjustReactionAndDiffusionRateForTemperature() {

void TsIRTConfiguration::AdjustReactionRateForPH(G4String pHOrConcentration) {
std::vector<G4double> AcidComponents;
std::vector<G4double> PhosphateComponents;
G4double Ionic = 0.0;
G4double HCon = 0.0;
G4double HCon25 = 1.00E-7;
G4double OHCon = 0.0;
G4double OHCon25 = 1.00E-7;
G4double HSO4Con = 0.0;
G4double PO4Con = 0.0;
G4double HPO4Con = 0.0;
G4double H2PO4Con = 0.0;
G4double H3PO4Con = 0.0;

if (pHOrConcentration == "PH" && fpHSolvent == "h2so4") {
AcidComponents = GetH2SO4ComponentsConcentrationPH(fpHValue);
Expand All @@ -1325,7 +1387,6 @@ void TsIRTConfiguration::AdjustReactionRateForPH(G4String pHOrConcentration) {
HSO4Con = AcidComponents[1];
OHCon = AcidComponents[3];
G4cout << "-- Adjust for PH of H2SO4 " << G4endl;

}

else if (pHOrConcentration == "Concentration" && fpHSolvent == "h2so4") {
Expand All @@ -1342,13 +1403,40 @@ void TsIRTConfiguration::AdjustReactionRateForPH(G4String pHOrConcentration) {
OHCon = 1E-14 / HCon;
AcidComponents = {HCon, 0.0, 0.0, OHCon, 0.0, 0.0};
Ionic = GetIonicStrength(AcidComponents);
G4cout << "-- Adjust for a generic substance " << G4endl;

} else {
G4cout << "-- Adjust for a generic substance with pH " << fpHValue << G4endl;
}

// PB: Assume standard values for phosphate buffer around pH7.4
else if (fpHSolvent == "PB" || fpHSolvent == "PBS" ) {
PhosphateComponents = GetPhosphateBufferComponentsConcentrationPHandIonicStrength(fpHValue);
Ionic = PhosphateComponents[0];
HCon = PhosphateComponents[1];
OHCon = PhosphateComponents[6];
PO4Con = PhosphateComponents[2];
HPO4Con = PhosphateComponents[3];
H2PO4Con = PhosphateComponents[4];
H3PO4Con = PhosphateComponents[5];
G4cout << "-- Adjust ionic strength for phosphate buffer at pH "<< fpHValue<< G4endl;

if (fpHSolvent == "PBS" ){
// Add monovalent salts from 1xPBS
// with physiological salts of 137 mM NaCl and 2.7mM KCl
G4double PBS_salts_conc = 139.7e-3;
Ionic = Ionic + GetIonicStrengthFromSalts(PBS_salts_conc,0.0,0.0);
G4cout << "-- Add 139.7mM monovalent salts for phosphate buffered saline (PBS)." << G4endl;
}
}

else {
G4cout << "-- Is doing nothing " << pHOrConcentration << " " << fpHSolvent << G4endl;
}

G4cout << G4endl;

G4cout << " ###--- Add Salt effects to Ionic strength ---###" << G4endl;
G4double additionalIonicFromSalts = GetIonicStrengthFromSalts(fMonovalentSalt,fDivalentSalt,fTrivalentSalt);
G4cout << " ###--- Add additional ionic strength from salts: " << additionalIonicFromSalts <<G4endl;
Ionic = Ionic + additionalIonicFromSalts;

G4cout << " ###-------- pH Scaling Starts ---------###" << G4endl;

for(size_t i = 0; i < fReactions.size(); i++) {
Expand Down Expand Up @@ -1465,6 +1553,16 @@ void TsIRTConfiguration::AdjustReactionRateForPH(G4String pHOrConcentration) {
setScavenging(fReactions[i].scavengingCapacity * OHCon / OHCon25);
} else if (ReactB == "HSO4^-1") {
setScavenging(fReactions[i].scavengingCapacity * HSO4Con);
} else if (ReactB == "H3PO4^0") {
setScavenging(fReactions[i].scavengingCapacity * H2PO4Con);
}
// Add dihydrogen phosphate Ion to the reactions with e_aq- + H2P04- 1.9 X 10^7/Ms (Buxton 343 at pH5-7)
else if (ReactB == "H2PO4^-1") {
setScavenging(fReactions[i].scavengingCapacity * H2PO4Con);
} else if (ReactB == "HPO4^-2") {
setScavenging(fReactions[i].scavengingCapacity * HPO4Con);
}else if (ReactB == "PO4^-3") {
setScavenging(fReactions[i].scavengingCapacity * PO4Con);
}
}

Expand Down
7 changes: 7 additions & 0 deletions processes/TsIRTConfiguration.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,10 @@ private:
G4String fpHSolvent;
G4double fpHSolventConcentration;
G4double fpHValue;
G4double fMonovalentSalt; // NaCl, KCl, ...
G4double fDivalentSalt; // MgCl2 ...
G4double fTrivalentSalt; // FeCl3 ...


G4bool fAllTotallyDiffusionControlled;

Expand Down Expand Up @@ -163,6 +167,9 @@ public:
std::vector<G4double> GetH2SO4ComponentsConcentrationP(G4double);
std::vector<G4double> GetH2SO4ComponentsConcentrationPH(G4double);
G4double GetIonicStrength(std::vector<G4double>);
G4double GetIonicStrengthFromSalts(G4double fMonovalentSalt, G4double fDivalentSalt, G4double fTrivalentSalt);
std::vector<G4double> GetPhosphateBufferComponentsConcentrationPHandIonicStrength(G4double pH);

inline std::map<G4String, G4int> GetMoleculeIDs() {return fMoleculesID;};
inline std::map<G4int, G4String> GetMoleculeNames() { return fMoleculesName;};
G4double IonicRate(G4double, TsMolecularReaction);
Expand Down
2 changes: 1 addition & 1 deletion userContributions/CellsWithNanoparticles/README.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ Examples can be found in the "examples/" directory.

4) Bugs:

Please report bugs to hahn@physik.fu-berlin.de or on https://github.com/BAMresearch/TOPAS-CellModels
Please report bugs to marc.hahn@uni-potsdam.de or on https://github.com/MarcBHahn/TOPAS-CellModels

5) Literature:

Expand Down
Loading