Add HIC Binding Models

.github/workflows/ci.yml

@@ -111,6 +111,7 @@ jobs:
           path: src
       - name: Install Dependencies
         run: |
+          sudo apt-get update
           sudo apt -y install \
             build-essential \
             libhdf5-dev \

doc/interface/binding/hic_constant_water_activity.rst

@@ -0,0 +1,78 @@
+.. _hic_constant_water_activity_config:
+
+HIC Constant Water Activity
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+**Group /input/model/unit_XXX/adsorption – ADSORPTION_MODEL = HIC_CONSTANT_WATER_ACTIVITY**
+
+
+``IS_KINETIC``
+   Selects kinetic or quasi-stationary adsorption mode: 1 = kinetic, 0 =
+   quasi-stationary. If a single value is given, the mode is set for all
+   bound states. Otherwise, the adsorption mode is set for each bound
+   state separately.
+
+===================  =========================  =======================
+**Type:** int        **Range:** {0,1}           **Length:** 1/NTOTALBND
+===================  =========================  =======================
+
+``HICCWA_KA``
+   Adsorption rate constant
+
+**Unit:** :math:`m_{MP}^{3}~m_{SP}^{-3}~s^{-1}`
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** NCOMP
+===================  =========================  =========================================
+
+``HICCWA_KD``
+   Desorption rate constant
+
+**Unit:** :math:`s^{-1}`
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** NCOMP
+===================  =========================  =========================================
+
+``HICCWA_NU``
+   Number of ligands per ligand-protein interaction
+
+**Unit: [-]**
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** NCOMP
+===================  =========================  =========================================
+
+
+``HICCWA_QMAX``
+   Maximum binding capacity
+
+**Unit:** :math:`mol~m_{SP}^{-3}`
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** NCOMP
+===================  =========================  =========================================
+
+
+``HICCWA_BETA0``
+   Parameters describing the number of highly ordered water molecules 
+   that stabilize the hydrophobic surfaces at infinitely diluted 
+   salt concentration
+
+**Unit: [-]**
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** 1
+===================  =========================  =========================================
+
+``HICCWA_BETA1``
+   Parameters describing the change in the number of highly ordered  
+   water molecules that stabilize the hydrophobic surfaces with
+   respect to changes in the salt concentration
+
+**Unit:** :math:`m_{MP}^{3}~mol^{-1}`
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** 1
+===================  =========================  =========================================
+

doc/interface/binding/hic_water_on_hydrophobic_surfaces.rst

@@ -0,0 +1,78 @@
+.. _hic_water_on_hydrophobic_surfaces_config:
+
+HIC Water on Hydrophobic Surfaces
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+**Group /input/model/unit_XXX/adsorption – ADSORPTION_MODEL = HIC_WATER_ON_HYDROPHOBIC_SURFACES**
+
+
+``IS_KINETIC``
+   Selects kinetic or quasi-stationary adsorption mode: 1 = kinetic, 0 =
+   quasi-stationary. If a single value is given, the mode is set for all
+   bound states. Otherwise, the adsorption mode is set for each bound
+   state separately.
+
+===================  =========================  =======================
+**Type:** int        **Range:** {0,1}           **Length:** 1/NTOTALBND
+===================  =========================  =======================
+
+``HICWHS_KA``
+   Adsorption rate constant
+
+**Unit:** :math:`m_{MP}^{3}~m_{SP}^{-3}~s^{-1}`
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** NCOMP
+===================  =========================  =========================================
+
+``HICWHS_KD``
+   Desorption rate constant
+
+**Unit:** :math:`s^{-1}`
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** NCOMP
+===================  =========================  =========================================
+
+``HICWHS_NU``
+   Number of ligands per ligand-protein interaction
+
+**Unit: [-]**
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** NCOMP
+===================  =========================  =========================================
+
+
+``HICWHS_QMAX``
+   Maximum binding capacity
+
+**Unit:** :math:`mol~m_{SP}^{-3}`
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** NCOMP
+===================  =========================  =========================================
+
+
+``HICWHS_BETA0``
+   Parameters describing the number of highly ordered water molecules 
+   that stabilize the hydrophobic surfaces at infinitely diluted 
+   salt concentration
+
+**Unit: [-]**
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** 1
+===================  =========================  =========================================
+
+``HICWHS_BETA1``
+   Parameters describing the change in the number of highly ordered  
+   water molecules that stabilize the hydrophobic surfaces with
+   respect to changes in the salt concentration
+
+**Unit:** :math:`m_{MP}^{3}~mol^{-1}`
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** 1
+===================  =========================  =========================================
+

doc/interface/binding/index.rst

@@ -65,4 +65,6 @@ This group also takes precedence over a possibly existing ``/input/model/unit_XX
     saska
     self_association
     kumar_langmuir
+    hic_water_on_hydrophobic_surfaces
+    hic_constant_water_activity
 

doc/literature.bib

@@ -413,4 +413,27 @@ @article{Benedikt2019
 doi = {10.1021/acs.est.8b05873},
 URL = { https://doi.org/10.1021/acs.est.8b05873},
 eprint = { https://doi.org/10.1021/acs.est.8b05873}
-}
+}
+@article{Wang2016,
+title = {Water on hydrophobic surfaces: Mechanistic modeling of hydrophobic interaction chromatography},
+journal = {Journal of Chromatography A},
+volume = {1465},
+pages = {71-78},
+year = {2016},
+issn = {0021-9673},
+doi = {https://doi.org/10.1016/j.chroma.2016.07.085},
+url = {https://www.sciencedirect.com/science/article/pii/S0021967316310263},
+author = {Gang Wang and Tobias Hahn and Jürgen Hubbuch},
+}
+@article{Jaepel2022,
+title = {Bayesian optimization using multiple directional objective functions allows the rapid inverse fitting of parameters for chromatography simulations},
+journal = {Journal of Chromatography A},
+volume = {1679},
+pages = {463408},
+year = {2022},
+issn = {0021-9673},
+doi = {https://doi.org/10.1016/j.chroma.2022.463408},
+url = {https://www.sciencedirect.com/science/article/pii/S0021967322005830},
+author = {Ronald Colin Jäpel and Johannes Felix Buyel},
+}
+

doc/modelling/binding/hic_constant_water_activity.rst

@@ -0,0 +1,18 @@
+.. _hic_constant_water_activity_model:
+
+HIC Constant Water Activity
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+This model implemments the HIC isotherm assuming a constant water activity as described by Jäpel and Buyel :cite:`Jaepel2022`.
+
+.. math::
+    \begin{align}
+        \beta &= \beta_0 e^{c_{p,0}\beta_1} \\
+        \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \left( 1 - \sum_j \frac{q_j}{q_{max,j}} \right)^{\nu_i} - k_{d,i} q_i 0.1^{\nu_i \beta}
+    \end{align}
+
+- Component :math:`c_0` is assumed to be salt without a bound state.
+- Multiple bound states are not supported.
+- Components without bound state (i.e., salt and non-binding components) are supported.
+
+For more information on model parameters required to define in CADET file format, see :ref:`hic_constant_water_activity_config`.
+

doc/modelling/binding/hic_water_on_hydrophobic_surfaces.rst

@@ -0,0 +1,20 @@
+.. _hic_water_on_hydrophobic_surfaces_model:
+
+HIC Water on Hydrophobic Surfaces
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This model implements a slightly modified version of the HIC isotherm by Wang et al. based on their 2016 paper :cite:`Wang2016`.
+A naive multicomponent version was added that reduces to the original formulation if only 1 binding species is present.
+
+.. math::
+
+    \begin{align}
+		\beta &= \beta_0 e^{c_{p,0}\beta_1} \\
+		\frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \left( 1 - \sum_j \frac{q_j}{q_{max,j}} \right)^{\nu_i} - k_{d,i} q_i  \left(\sum_j q_j \right)^{\nu_i \beta}
+    \end{align}
+
+- Component :math:`c_0` is assumed to be salt without a bound state.
+- Multiple bound states are not supported.
+- Components without bound state (i.e., salt and non-binding components) are supported.
+
+For more information on model parameters required to define in CADET file format, see :ref:`hic_water_on_hydrophobic_surfaces_config`.

doc/modelling/binding/index.rst

@@ -257,6 +257,16 @@ The models also differ in whether a mobile phase modifier (e.g., salt) is suppor
      - ✓
      - ✓
      - ✓
+   * - :ref:`hic_water_on_hydrophobic_surfaces_model`
+     - ✓
+     - x
+     - ✓
+     - x
+   * - :ref:`hic_constant_water_activity_model`
+     - ✓
+     - x
+     - ✓
+     - x
 
 
 .. toctree::

src/libcadet/BindingModelFactory.cpp

@@ -42,6 +42,8 @@ namespace cadet
 			void registerLangmuirLDFModel(std::unordered_map<std::string, std::function<model::IBindingModel* ()>>& bindings);
 			void registerLangmuirLDFCModel(std::unordered_map<std::string, std::function<model::IBindingModel* ()>>& bindings);
 			void registerBiLangmuirLDFModel(std::unordered_map<std::string, std::function<model::IBindingModel* ()>>& bindings);
+			void registerHICWaterOnHydrophobicSurfacesModel(std::unordered_map<std::string, std::function<model::IBindingModel*()>>& bindings);
+			void registerHICConstantWaterActivityModel(std::unordered_map<std::string, std::function<model::IBindingModel*()>>& bindings);
 		}
 	}
 
@@ -69,6 +71,8 @@ namespace cadet
 		model::binding::registerLangmuirLDFModel(_bindingModels);
 		model::binding::registerLangmuirLDFCModel(_bindingModels);
 		model::binding::registerBiLangmuirLDFModel(_bindingModels);
+		model::binding::registerHICWaterOnHydrophobicSurfacesModel(_bindingModels);
+		model::binding::registerHICConstantWaterActivityModel(_bindingModels);
 		registerModel<model::SimplifiedMultiStateStericMassActionBinding>();
 	}
 

src/libcadet/CMakeLists.txt

@@ -98,6 +98,8 @@ set(LIBCADET_BINDINGMODEL_SOURCES
 	${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/LangmuirLDFBinding.cpp
 	${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/LangmuirLDFCBinding.cpp
 	${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/BiLangmuirLDFBinding.cpp
+	${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICWaterOnHydrophobicSurfacesBinding.cpp
+	${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICConstantWaterActivityBinding.cpp
 )
 
 # LIBCADET_REACTIONMODEL_SOURCES holds all source files of reaction models

src/libcadet/model/Parameters.hpp

@@ -132,6 +132,24 @@ class ScalarParameter
 		*_p = paramProvider.getDouble(varName);
 	}
 
+	/**
+	 * @brief Reads parameters and verifies them
+	 * @details See IBindingModel::configure() for details.
+	 * @param [in] varName Name of the parameter
+	 * @param [in] paramProvider IParameterProvider used for reading parameters
+	 * @param [in] nComp Number of components
+	 * @param [in] nBoundStates Array with number of bound states for each component
+	 * @param [in] defaultVal Default value if parameter is not available in @p paramProvider
+	 * @return @c true if the parameters were read and validated successfully, otherwise @c false
+	 */
+	inline void configure(const std::string& varName, IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates, double defaultVal)
+	{
+		if (paramProvider.exists(varName))
+			*_p = paramProvider.getDouble(varName);
+		else
+			*_p = defaultVal;
+	}
+
 	/**
 	 * @brief Registers the parameters in a map for further use
 	 * @param [in] varName Name of the parameter

src/libcadet/model/binding/ExternalFunctionTemplate.cpp

@@ -94,9 +94,17 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase
 	{% for p in constantParameters %}
 		{% if not existsIn(p, "skipConfig") %}
 			{% if length(p/varName) > 1 %}
-				_{{ p/objName }}.configure("{{ p/confPrefix }}", paramProvider, nComp, nBoundStates);
+				{% if existsIn(p, "default") %}
+					_{{ p/objName }}.configure("{{ p/confPrefix }}", paramProvider, nComp, nBoundStates, {{ p/default }});
+				{% else %}
+					_{{ p/objName }}.configure("{{ p/confPrefix }}", paramProvider, nComp, nBoundStates);
+				{% endif %}
 			{% else %}
-				_{{ p/varName }}.configure("{{ p/confName }}", paramProvider, nComp, nBoundStates);
+				{% if existsIn(p, "default") %}
+					_{{ p/varName }}.configure("{{ p/confName }}", paramProvider, nComp, nBoundStates, {{ p/default }});
+				{% else %}
+					_{{ p/varName }}.configure("{{ p/confName }}", paramProvider, nComp, nBoundStates);
+				{% endif %}
 			{% endif %}
 		{% endif %}
 	{% endfor %}
@@ -265,9 +273,17 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase
 	{% for p in constantParameters %}
 		{% if not existsIn(p, "skipConfig") %}
 			{% if length(p/varName) > 1 %}
-				_{{ p/objName }}.configure("{{ p/confPrefix }}", paramProvider, nComp, nBoundStates);
+				{% if existsIn(p, "default") %}
+					_{{ p/objName }}.configure("{{ p/confPrefix }}", paramProvider, nComp, nBoundStates, {{ p/default }});
+				{% else %}
+					_{{ p/objName }}.configure("{{ p/confPrefix }}", paramProvider, nComp, nBoundStates);
+				{% endif %}
 			{% else %}
-				_{{ p/varName }}.configure("{{ p/confName }}", paramProvider, nComp, nBoundStates);
+				{% if existsIn(p, "default") %}
+					_{{ p/varName }}.configure("{{ p/confName }}", paramProvider, nComp, nBoundStates, {{ p/default }});
+				{% else %}
+					_{{ p/varName }}.configure("{{ p/confName }}", paramProvider, nComp, nBoundStates);
+				{% endif %}
 			{% endif %}
 		{% endif %}
 	{% endfor %}