Skip to content

Commit c6880be

Browse files
ronald-jaepelronald-jaepel
committed
Add pH dependence
1 parent 96f2bea commit c6880be

File tree

1 file changed

+54
-20
lines changed

1 file changed

+54
-20
lines changed

‎src/libcadet/model/binding/HICUnifiedBinding.cpp

Lines changed: 54 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -35,15 +35,18 @@
3535
{ "type": "ScalarParameter", "varName": "beta0", "confName": "HICUNI_BETA0"},
3636
{ "type": "ScalarParameter", "varName": "beta1", "confName": "HICUNI_BETA1"},
3737
{ "type": "ScalarComponentDependentParameter", "varName": "kA", "confName": "HICUNI_KA"},
38+
{ "type": "ScalarComponentDependentParameter", "varName": "kALin", "confName": "HICUNI_KA_LIN"},
3839
{ "type": "ScalarComponentDependentParameter", "varName": "kD", "confName": "HICUNI_KD"},
3940
{ "type": "ScalarComponentDependentParameter", "varName": "kP", "confName": "HICUNI_KP"},
4041
{ "type": "ScalarComponentDependentParameter", "varName": "kS", "confName": "HICUNI_KS"},
4142
{ "type": "ScalarComponentDependentParameter", "varName": "epsilon", "confName": "HICUNI_EPSILON"},
4243
{ "type": "ScalarComponentDependentParameter", "varName": "nu", "confName": "HICUNI_NU"},
44+
{ "type": "ScalarComponentDependentParameter", "varName": "nuLin", "confName": "HICUNI_NU_LIN"},
4345
{ "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "HICUNI_QMAX"}
4446
],
4547
"constantParameters":
4648
[
49+
{ "type": "ScalarParameter", "varName": "refPh", "default": 7.0, "objName": "referencePh", "confName": "HICUNI_PH_REF"},
4750
{ "type": "ScalarParameter", "varName": "rho", "default": 0.00003348, "objName": "osmoticEffect", "confName": "HICUNI_RHO"}
4851
4952
]
@@ -55,13 +58,16 @@
5558
beta0 = bulk-like ordered water molecules at infinitely diluted salt concentration
5659
beta1 = influence of salt concentration on bulk-like ordered water molecules
5760
kA = Adsorption constant
61+
kALin = Linear dependency of ka on pH
5862
kD = Desorption constant
5963
kP = Influence of protein concentration on protein activity
6064
kS = Influence of salt concentration on protein activity
6165
epsilon = Influence of bound protein concentration on bound protein activity
6266
nu = Number of binding sites
67+
nuLin = Linear dependency of nu on pH
6368
qMax = Maximum binding capacity
6469
rh0 = osmotic effect of salt concentration on water activity. Calculated as osmotic_coefficient * molar_weight_of_water * ion_number in
70+
refpH = reference pH (defaults to 7)
6571
*/
6672

6773
namespace cadet
@@ -75,8 +81,10 @@ namespace cadet
7581
inline bool HICUNIFIEDParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
7682
{
7783
if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size())
78-
|| (_kA.size() != _kP.size()) || (_kA.size() != _kS.size()) || (_kA.size() != _epsilon.size()) || (_kA.size() < nComp))
79-
throw InvalidParameterException("HICUNI_KA, HICUNI_KD, HICUNI_NU, and HICUNI_QMAX have to have the same size");
84+
|| (_kA.size() != _kP.size()) || (_kA.size() != _kS.size()) || (_kA.size() != _epsilon.size())
85+
|| (_kA.size() != _nuLin.size()) || (_kA.size() != _kALin.size())
86+
|| (_kA.size() < nComp))
87+
throw InvalidParameterException("HICUNI_KA, HICUNI_KD, HICUNI_NU, and HICUNI_QMAX, HICUNI_KA_LIN, HICUNI_NU_LIN have to have the same size");
8088

8189
return true;
8290
}
@@ -86,8 +94,10 @@ namespace cadet
8694
inline bool ExtHICUNIFIEDParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
8795
{
8896
if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size())
89-
|| (_kA.size() != _kP.size()) || (_kA.size() != _kS.size()) || (_kA.size() != _epsilon.size()) || (_kA.size() < nComp))
90-
throw InvalidParameterException("KA, KD, NU, and QMAX have to have the same size");
97+
|| (_kA.size() != _kP.size()) || (_kA.size() != _kS.size()) || (_kA.size() != _epsilon.size())
98+
|| (_kA.size() != _nuLin.size()) || (_kA.size() != _kALin.size())
99+
|| (_kA.size() < nComp))
100+
throw InvalidParameterException("HICUNI_KA, HICUNI_KD, HICUNI_NU, and HICUNI_QMAX, HICUNI_KA_LIN, HICUNI_NU_LIN have to have the same size");
91101

92102
return true;
93103
}
@@ -97,7 +107,10 @@ namespace cadet
97107
* @brief Defines a unified HIC isotherm combining publications by Deichert et al 2010, Mollerup et al 2006, and Jäpel and Buyel 2022
98108
* @details Implements a unified HIC isotherm : \f[ \begin{align}
99109
* \beta &= \beta_0 e^{c_{p,0}\beta_1}\\
100-
* \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i}e^{k_{p,i} c_{p,i}+k_{s,i} c_{s}}(1-\sum{\frac{q_p}{q_{max}}})^{ν_i} -k_{d,i} q_{p,i}(1+\epsilon q_{p,i}) (e^{ρc_{s}})^{ν_i β_0 e^{β_1 c_s}}
110+
* k_{a,i} &= k_{a,i,0} \exp\left( k_{a,i,lin} (\mathrm{pH}-\mathrm{pH}_{ref})\right)\\
111+
* \nu_i &= \nu_{i,0} + \nu_{i,lin} (\mathrm{pH}-\mathrm{pH}_{ref})\\
112+
* \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i}e^{k_{p,i} c_{p,i}+k_{s,i} c_{s}}(1-\sum{\frac{q_p}{q_{max}}})^{\nu_i}
113+
* -k_{d,i} q_{p,i}(1+\epsilon q_{p,i}) (e^{\rho c_{s}})^{\nu_i \beta_0 e^{\beta_1 c_s}}
101114
* \end{align} \f]
102115
* Component @c 0 is assumed to be salt without a bound state. Multiple bound states are not supported.
103116
* Components without bound state (i.e., salt and non-binding components) are supported.
@@ -122,6 +135,10 @@ namespace cadet
122135
if (nBound[0] != 0)
123136
throw InvalidParameterException("HICUNIFIED binding model requires exactly zero bound states for salt component");
124137

138+
139+
// Guarantee that modifier component is non-binding
140+
if (nBound[1] != 0)
141+
throw InvalidParameterException("HICUNIFIED binding model requires non-binding pH modifier component (NBOUND[1] = 0)");
125142
return res;
126143
}
127144

@@ -167,6 +184,10 @@ namespace cadet
167184

168185
typename ParamHandler_t::ParamsHandle const _p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace);
169186

187+
// Pseudo component 1 is pH
188+
const CpStateType pH = yCp[1];
189+
const ParamType refPh = static_cast<ParamType>(_p->refPh);
190+
170191
const ParamType beta0 = static_cast<ParamType>(_p->beta0);
171192
const ParamType beta1 = static_cast<ParamType>(_p->beta1);
172193

@@ -183,7 +204,7 @@ namespace cadet
183204
StateParamType qSum = 1.0;
184205

185206
unsigned int bndIdx = 0;
186-
for (int i = 0; i < _nComp; ++i)
207+
for (int i = 2; i < _nComp; ++i)
187208
{
188209
// Skip components without bound states (bound state index bndIdx is not advanced)
189210
if (_nBoundStates[i] == 0)
@@ -196,23 +217,26 @@ namespace cadet
196217
}
197218

198219
bndIdx = 0;
199-
for (int i = 0; i < _nComp; ++i)
220+
for (int i = 2; i < _nComp; ++i)
200221
{
201222
// Skip components without bound states (bound state index bndIdx is not advanced)
202223
if (_nBoundStates[i] == 0)
203224
continue;
204225

205226
const ParamType kD = static_cast<ParamType>(_p->kD[i]);
206227
const ParamType kA = static_cast<ParamType>(_p->kA[i]);
228+
const ParamType kALin = static_cast<ParamType>(_p->kALin[i]);
207229
const ParamType kP = static_cast<ParamType>(_p->kP[i]);
208230
const ParamType kS = static_cast<ParamType>(_p->kS[i]);
209231
const ParamType nu = static_cast<ParamType>(_p->nu[i]);
232+
const ParamType nuLin = static_cast<ParamType>(_p->nuLin[i]);
210233
const ParamType epsilon = static_cast<ParamType>(_p->epsilon[i]);
211234

212-
/*res[bndIdx] = kD * y[bndIdx] * pow(water_activity, nu * beta)
213-
- kA * pow(qSum, nu) * yCp[i];*/
214-
res[bndIdx] = kD * y[bndIdx] * (1 + epsilon * y[bndIdx]) * pow(water_activity, nu * beta)
215-
- kA * pow(qSum, nu) * yCp[i] * exp(kP * yCp[i] + kS * yCp[0]);
235+
const CpStateParamType kApH = kA * exp(kALin * (pH - refPh));
236+
const CpStateParamType nupH = nu + nuLin * (pH - refPh);
237+
238+
res[bndIdx] = kD * y[bndIdx] * (1 + epsilon * y[bndIdx]) * pow(water_activity, nupH * beta)
239+
- kApH * pow(qSum, nupH) * yCp[i] * exp(kP * yCp[i] + kS * yCp[0]);
216240

217241
// Next bound component
218242
++bndIdx;
@@ -228,6 +252,10 @@ namespace cadet
228252

229253
using std::pow, std::exp;
230254

255+
// Pseudo component 1 is pH
256+
const double pH = yCp[1];
257+
const double refPh = static_cast<double>(_p->refPh);
258+
231259
const double beta0 = static_cast<double>(_p->beta0);
232260
const double beta1 = static_cast<double>(_p->beta1);
233261

@@ -239,7 +267,7 @@ namespace cadet
239267
double qSum = 1.0;
240268

241269
unsigned int bndIdx = 0;
242-
for (int i = 0; i < _nComp; ++i)
270+
for (int i = 2; i < _nComp; ++i)
243271
{
244272
// Skip components without bound states (bound state index bndIdx is not advanced)
245273
if (_nBoundStates[i] == 0)
@@ -252,7 +280,7 @@ namespace cadet
252280
}
253281

254282
bndIdx = 0;
255-
for (int i = 0; i < _nComp; ++i)
283+
for (int i = 2; i < _nComp; ++i)
256284
{
257285
// Skip components without bound states (bound state index bndIdx is not advanced)
258286
if (_nBoundStates[i] == 0)
@@ -264,15 +292,21 @@ namespace cadet
264292
const double kS = static_cast<double>(_p->kS[i]);
265293
const double nu = static_cast<double>(_p->nu[i]);
266294
const double epsilon = static_cast<double>(_p->epsilon[i]);
295+
const double kALin = static_cast<double>(_p->kALin[i]);
296+
const double nuLin = static_cast<double>(_p->nuLin[i]);
297+
298+
const double kApH = kA * exp(kALin * (pH - refPh));
299+
const double nupH = nu + nuLin * (pH - refPh);
300+
267301

268302
// dres_i / dc_{p,0}
269-
jac[-bndIdx - offsetCp] = kD * y[bndIdx] * (1 + epsilon * y[bndIdx]) * (-rho) * nu * beta0 * (beta1 * yCp[0] + 1) *
270-
exp(beta1 * yCp[0] - rho * nu * beta0 * exp(beta1 * yCp[0]) * yCp[0])
271-
- kA * yCp[i] * pow(qSum, nu) * exp(kP * yCp[i] + kS * yCp[0]) * kS;
303+
jac[-bndIdx - offsetCp] = kD * y[bndIdx] * (1 + epsilon * y[bndIdx]) * (-rho) * nupH * beta0 * (beta1 * yCp[0] + 1) *
304+
exp(beta1 * yCp[0] - rho * nupH * beta0 * exp(beta1 * yCp[0]) * yCp[0])
305+
- kApH * yCp[i] * pow(qSum, nupH) * exp(kP * yCp[i] + kS * yCp[0]) * kS;
272306

273307
// dres_i / dc_{p,i}
274-
jac[i - bndIdx - offsetCp] = -kA * pow(qSum, nu) * exp(kP * yCp[i] + kS * yCp[0])
275-
- kA * yCp[i] * pow(qSum, nu) * exp(kP * yCp[i] + kS * yCp[0]) * kP;
308+
jac[i - bndIdx - offsetCp] = -kApH * pow(qSum, nupH) * exp(kP * yCp[i] + kS * yCp[0])
309+
- kApH * yCp[i] * pow(qSum, nupH) * exp(kP * yCp[i] + kS * yCp[0]) * kP;
276310
// Getting to c_{p,i}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +i to c_{p,i}.
277311
// This means jac[i - bndIdx - offsetCp] corresponds to c_{p,i}.
278312

@@ -286,14 +320,14 @@ namespace cadet
286320
continue;
287321

288322
// dres_i / dq_j
289-
jac[bndIdx2 - bndIdx] = kA * yCp[i] * nu * pow(qSum, nu - 1) / (static_cast<double>(_p->qMax[j])) * exp(kP * yCp[i] + kS * yCp[0]);
323+
jac[bndIdx2 - bndIdx] = kApH * yCp[i] * nupH * pow(qSum, nupH - 1) / (static_cast<double>(_p->qMax[j])) * exp(kP * yCp[i] + kS * yCp[0]);
290324
// Getting to q_j: -bndIdx takes us to q_0, another +bndIdx2 to q_j. This means jac[bndIdx2 - bndIdx] corresponds to q_j.
291325

292326
++bndIdx2;
293327
}
294328

295329
// Add to dres_i / dq_i
296-
jac[0] += (1 + 2 * epsilon * y[bndIdx]) * kD * pow(water_activity, nu * beta);
330+
jac[0] += (1 + 2 * epsilon * y[bndIdx]) * kD * pow(water_activity, nupH * beta);
297331

298332
// Advance to next flux and Jacobian row
299333
++bndIdx;

0 commit comments

Comments
 (0)