CADET isotherm model

I was exploring CADET different isotherm models. Currently, there is no HIC model or any pH-related model on the list here. Am I right? Binding models — CADET

However, there is an example in the CADET-Python model named as HICWang. I assume this is a HIC model associated with his paper in 2016. The code doesn’t work well on my end (SMAmodel does work) . I don’t know if the adsorption model is not correct since HICWANG is not listed in the binding model list or if anything else is broken.

root.input.model.unit_001.adsorption_model = ‘HICWANG’’

Is there a way that we can include more isotherm models or work together on that? I do recall seeing publications using CADET with a different model not listed in the summary table.

Thanks ahead for your help.

Best,
Yiran

Hi Yiran,

to my knowledge there is currently no HIC model. I assume the example in the CADET-Python package was some custom model. Maybe @w.heymann can comment on this?

In case you’re interested, the implementation of binding models in the source code lives here.
We’re always happy to assist if you’re interested in implementing a new model. In fact, @j.hassan is currently working on a “How-to” guide for this exact purpose. Maybe you’d be willing to test the guide and provide some feedback?

Best

Johannes

Hey,

I’ve put the HICWANG code (written by @w.heymann & edited by me) into a pull request on github.

Best wishes,
Ron

1 Like

Also, I think I should add that the reason the Wang Isotherm wasn’t added yet is that it shows unrealistic behavior. We have two alternative versions queued for publication and I can update you once they are available.

1 Like

Not only was it unrealistic the units on the equation didn’t even work out. :slight_smile:

Hi Ron, just want to follow up about the HIC isotherm model :slight_smile:

There are, as of this morning, two HIC models in the main branch of CADET. Do you know how to compile CADET from source?

The first isotherm is the one published by Gang Wang et al. in “Water on Hydrophobic Surfaces” (2016) and is called HIC_WHS based on the paper title. As mentioned above, this one exhibits a few unexpected behaviors.

The second is a simplification of the WHS isotherm that solves the unexpected behaviors and is called HIC_CWA for “Constant Water Activity”.

1 Like

It might be worth looking into the HIC model used by J.M. Mollerup, T.B. Hansen, S. Kidal, A. Staby: Quality by design-Thermodynamic modelling of chromatographic separation of proteins. In this paper, the isotherm is only written in its equilibrium form, so I have expressed this in the kinetic form shown below

\frac{\partial q_i}{\partial t}=k_{a, i} \left(1-\sum_{j=1}^N \frac{q_j}{q_{\max , j}}\right)^{n_i} \tilde{\gamma}_i c_{p, i}-k_{d, i}q_i

This has been employed in several publications and includes the asymmetric activity coefficient

\tilde{\gamma}_i=\frac{\gamma_i}{\gamma^{\infty, \omega}} \approx \exp \left(K_{p, i} c_{p, i}+K_{s, i} c_s\right)

which is a very powerful term—see our recent paper in JChromA on isotherm models in multimodal chromatography, not HIC but I think the takeaways on modeling hydrophobic interactions are relevant here—and does not include the \beta_i, which is redundant with K_{s,i} in its dependence on salt concentration c_s and comparatively a less useful term.

The one thing I don’t like about this formulation is that it does not include ligand density, but I think for now it is sufficient—especially because ligand density for HIC resins is very rarely measured or provided by the manufacturer. If this is included, it should come along with the steric factor for hydrophobic interactions s_i, which could be pH-dependent. Lastly, pH could be included here by adding the dependence for k_{a,i} and/or k_{d,i}

k_{a,i} = k_{a,i,0} \exp \left(k_{a,i,1}\mathrm{pH} + k_{a,i,2}\mathrm{pH}^2 \right)

to the adsorption term in the \frac{\partial q_i}{\partial t} formulation. The asymmetric activity coefficient has been shown by Mollerup to be pH-dependent, but this need not be included as to avoid redundancy with the aforementioned pH-dependent affinity term.

The pH-dependent parameters and activity coefficient terms can be included in the k_{d,i} term for generality, but I think actually using these would most certainly lead to overparameterization of the isotherm. Actually, I always set k_{d,i} to unity such that k_{a,i} is equivalent to the equilibrium constant.

The stoichiometric term for hydrophobic interactions n_i can also be pH dependent and can be expressed analogously to \nu_i in the GIEX model.

n_i = n_{i,0} + n_{i,1}\mathrm{pH} + n_{i,2}\mathrm{pH}^2

If I could make a feature request, I think this would be a useful one. This model can be used for HIC or for multimodal in the pure HIC regime—when binding affinity and capacity are strictly increasing (electrostatic interactions are negligible) with respect to salt concentration, such as multimodal anion exchange for mAbs with pH << pI. Otherwise, the GIEX model can be employed for multimodal and, of course, for IEC. This way, we have two isotherms in CADET that are powerful and particularly general in their scope. We could call this one the generalized hydrophobic interaction (GHI) isotherm or something like that.

Thoughts?

Thank you so much for the update! I would love to try it out!

Compiling from source is kind of tricky. I have tried in the past following the instruction without success.

Alright, which platform (Windows/Linux/etc.) do you need binaries for?

I want to put a paper on HIC isotherms together and will see if I can add this isotherm to the mix.

2 Likes

I use Windows. Thanks!!!

CADET-commit 4eeeed5-win-x64.zip (3.9 MB)
Let me know if you run into any problems.

1 Like

I got this error message.

conda install -c conda-forge hdf5=1.12.2 doesn’t work. It said it conflicted with multiple current packages.

CompletedProcess(args=[‘C:\Users\wangy282\OneDrive - Bristol Myers Squibb\Documents\CADETYW\CADET-commit 4eeeed5-win-x64\bin\cadet-cli’, ‘model.h5’], returncode=3221226505, stdout=b’', stderr=b"Warning! HDF5 library version mismatched error\r\nThe HDF5 header files used to compile this application do not match\r\nthe version used by the HDF5 library to which this application is linked.\r\nData corruption or segmentation faults may occur if the application continues.\r\nThis can happen when an application was compiled by one version of HDF5 but\r\nlinked with a different version of static or shared HDF5 library.\r\nYou should recompile the application or check your shared library related\r\nsettings such as ‘LD_LIBRARY_PATH’.\r\nYou can, at your own risk, disable this warning by setting the environment\r\nvariable ‘HDF5_DISABLE_VERSION_CHECK’ to a value of ‘1’.\r\nSetting it to 2 or higher will suppress the warning messages totally.\r\nHeaders are 1.12.2, library is 1.10.6\r\n SUMMARY OF THE HDF5 CONFIGURATION\r\n =================================\r\n\r\nGeneral Information:\r\n-------------------\r\n HDF5 Version: 1.10.6\r\n Configured on: 2022-06-14\r\n Configured by: NMake Makefiles\r\n Host system: Windows-10.0.17763\r\n Uname information: Windows\r\n Byte sex: little-endian\r\n Installation point: C:/ci/hdf5_1655187604574/_h_env/Library\r\n\r\nCompiling Options:\r\n------------------\r\n Build Mode: RELEASE\r\n Debugging Symbols: OFF\r\n Asserts: OFF\r\n Profiling: OFF\r\n Optimization Level: OFF\r\n\r\nLinking Options:\r\n----------------\r\n Libraries: \r\n Statically Linked Executables: OFF\r\n LDFLAGS: /machine:x64\r\n H5_LDFLAGS: \r\n AM_LDFLAGS: \r\n Extra libraries: \r\n Archiver: C:/Program Files (x86)/Microsoft Visual Studio/2017/Professional/VC/Tools/MSVC/14.16.27023/bin/Hostx64/x64/lib.exe\r\n Ranlib: :\r\n\r\nLanguages:\r\n----------\r\n C: YES\r\n C Compiler: C:/Program Files (x86)/Microsoft Visual Studio/2017/Professional/VC/Tools/MSVC/14.16.27023/bin/Hostx64/x64/cl.exe 19.16.27043.0\r\n CPPFLAGS: \r\n H5_CPPFLAGS: \r\n AM_CPPFLAGS: \r\n CFLAGS: /DWIN32 /D_WINDOWS /W3\r\n H5_CFLAGS: \r\n AM_CFLAGS: \r\n Shared C Library: YES\r\n Static C Library: YES\r\n\r\n Fortran: OFF\r\n Fortran Compiler: \r\n Fortran Flags: \r\n H5 Fortran Flags: \r\n AM Fortran Flags: \r\n Shared Fortran Library: YES\r\n Static Fortran Library: YES\r\n\r\n C++: ON\r\n C++ Compiler: C:/Program Files (x86)/Microsoft Visual Studio/2017/Professional/VC/Tools/MSVC/14.16.27023/bin/Hostx64/x64/cl.exe 19.16.27043.0\r\n C++ Flags: /DWIN32 /D_WINDOWS /W3 /GR /EHsc\r\n H5 C++ Flags: \r\n AM C++ Flags: \r\n Shared C++ Library: YES\r\n Static C++ Library: YES\r\n\r\n JAVA: OFF\r\n JAVA Compiler: \r\n\r\nFeatures:\r\n---------\r\n Parallel HDF5: OFF\r\nParallel Filtered Dataset Writes: \r\n Large Parallel I/O: \r\n High-level library: ON\r\n Build HDF5 Tests: ON\r\n Build HDF5 Tools: ON\r\n Threadsafety: ON\r\n Default API mapping: v110\r\n With deprecated public symbols: ON\r\n I/O filters (external): DEFLATE\r\n MPE: \r\n Direct VFD: \r\n (Read-Only) S3 VFD: \r\n (Read-Only) HDFS VFD: \r\n dmalloc: \r\n Packages w/ extra debug output: \r\n API Tracing: OFF\r\n Using memory checker: OFF\r\n Memory allocation sanity checks: OFF\r\n Function Stack Tracing: OFF\r\n Strict File Format Checks: OFF\r\n Optimization Instrumentation: \r\nBye…\r\n")
Traceback (most recent call last):

File “C:\Users\wangy282\OneDrive - Bristol Myers Squibb\Documents\CADETYW\Python Example\230505_HIC_NFPB.py”, line 146, in
raise Exception(“Simulation failed”)

Exception: Simulation failed

Hi,

have you tried a new, clean conda environment?

best

Johannes

Try this: CADET-commit 4eeeed5-win-x64.zip (5.2 MB) It includes a matching hdf5.dll

1 Like

You rock! I am able to get the constant water model coverage. Maybe need to tweak the parameters but water on hydrophobic surface is always fail.

I do observe the model is not responsive to the change of modifier. The manual (Binding models — CADET) also indicates that model is not responsive to salt. But there is cp,0 in the equation. Am I miss something here?

1 Like

Which model parameters are you using?

And one thing to note is, that CADET assumes the salt to be the first component, so if you have your protein as the first model parameter and the salt as the second, it won’t work.

I am testing the parameters from Wang et al 2016.

model = Cadet()
model.root.input.model.nunits = 3
###Inlet
model.root.input.model.unit_000.unit_type = 'INLET'
model.root.input.model.unit_000.ncomp = 2
model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'

#ColumnLoad
model.root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL'
model.root.input.model.unit_001.ncomp = 2

## Geometry
model.root.input.model.unit_001.col_length = 0.059                # m
model.root.input.model.unit_001.cross_section_area = 0.000034       # m
model.root.input.model.unit_001.col_porosity = 0.35              # -
model.root.input.model.unit_001.par_porosity =0.91               # -
model.root.input.model.unit_001.par_radius = 2.1e-5              # m

                                                              
## Transport
model.root.input.model.unit_001.col_dispersion = 5.6e-8          # m^2 / s (interstitial volume)
model.root.input.model.unit_001.film_diffusion = [1.0E-6,1.0E-6]         # m / s
model.root.input.model.unit_001.par_diffusion = [3e-9,3e-11]        # m^2 / s (mobile phase)  
model.root.input.model.unit_001.par_surfdiffusion = [0.0,0]      # m^2 / s (solid phase)



# HIC Binding
model.root.input.model.unit_001.adsorption_model = 'HIC_CONSTANT_WATER_ACTIVITY'
model.root.input.model.unit_001.adsorption.is_kinetic = False    # Kinetic binding
model.root.input.model.unit_001.adsorption.HICCWA_KA = [0, 6e4]      # m^3 / (mol * s)   (mobile phase)
model.root.input.model.unit_001.adsorption.HICCWA_KD = [0, 2e2]      # 1 / s (desorption)
model.root.input.model.unit_001.adsorption.HICCWA_NU = [0,1]  
model.root.input.model.unit_001.adsorption.HICCWA_QMAX = [0, 4]  # mol / m^3 
model.root.input.model.unit_001.adsorption.HICCWA_BETA0 = [0,100] 
model.root.input.model.unit_001.adsorption.HICCWA_BETA1 = [0,10]  
model.root.input.model.unit_001.init_c = [000.0,0] #low salt
#model.root.input.model.unit_001.init_c = [1000.0,0] #high salt
model.root.input.model.unit_001.init_q = [0.0,0] 

### Grid cells
model.root.input.model.unit_001.discretization.ncol = 20
model.root.input.model.unit_001.discretization.npar = 5

### Bound states
model.root.input.model.unit_001.discretization.nbound = [0,1]

### Other options
model.root.input.model.unit_001.discretization.par_disc_type = 'EQUIDISTANT_PAR'    
model.root.input.model.unit_001.discretization.use_analytic_jacobian = 1
model.root.input.model.unit_001.discretization.reconstruction = 'WENO'
model.root.input.model.unit_001.discretization.gs_type = 1
model.root.input.model.unit_001.discretization.max_krylov = 0
model.root.input.model.unit_001.discretization.max_restarts = 10
model.root.input.model.unit_001.discretization.schur_safety = 1.0e-8

model.root.input.model.unit_001.discretization.weno.boundary_model = 0
model.root.input.model.unit_001.discretization.weno.weno_eps = 1e-10
model.root.input.model.unit_001.discretization.weno.weno_order = 3

###ColumnOutlet
model.root.input.model.unit_002.unit_type = 'OUTLET'
model.root.input.model.unit_002.ncomp = 2

### Section
model.root.input.solver.sections.nsec = 2
model.root.input.solver.sections.section_times = [0.0, 20.5/.26*60,40/.26*60]   # s
model.root.input.solver.sections.section_continuity = [0,]

#Load
model.root.input.model.unit_000.sec_000.const_coeff = [000,13.18/150] # mol / m^3 #lowsalt
model.root.input.model.unit_000.sec_000.const_coeff = [1000,13.18/150] # mol / m^3 #highsalt
model.root.input.model.unit_000.sec_000.lin_coeff = [0.0,0]
model.root.input.model.unit_000.sec_000.quad_coeff = [0.0,0]
model.root.input.model.unit_000.sec_000.cube_coeff = [0.0,0]

#Chase
model.root.input.model.unit_000.sec_001.const_coeff = [0,0] # mol / m^3

model.root.input.model.connections.nswitches = 1
model.root.input.model.connections.switch_000.section = 0
model.root.input.model.connections.switch_000.connections = [
    0, 1, -1, -1, 4.3/1e9,  # [unit_000, unit_001, all components, all components, Q/ m^3*s^-1 
    1, 2, -1, -1, 4.3/1e9]  # [unit_001, unit_002, all components, all components, Q/ m^3*s^-1 
model.root.input.model.solver.gs_type = 1
model.root.input.model.solver.max_krylov = 0
model.root.input.model.solver.max_restarts = 10
model.root.input.model.solver.schur_safety = 1e-8

# Number of cores for parallel simulation
model.root.input.solver.nthreads = 3

Alright, first thing to note is that parameters from Wang for the WHS aren’t going to work well with the CWA isotherm.

Secondly, Wang et al use mol/L while CADET uses the SI units of mol/m^3 (i.e. mmol/L). Therefore all parameters that get multiplied with a concentration need to be divided by 1000 going from Wang to CADET. A \beta_1 of 10 multiplied with 1000 mmol/L salt is going to get way too big and lead to problems. I’ve added param values that work for us below, but those are estimated with a salt concentration going from 3000 mmol/L to 0 mmol/L, not 1000 to 0 like in your example.

Lastly, the dimensionality of \beta_0 and \beta_1 should be 1 value total, instead of one value per component. That is probably why you didn’t get any salt influence, because the \beta_1 that CADET read was the 0 in first position, so 0 * c_s.

HIC Binding CWA

model.root.input.model.unit_001.adsorption_model = ‘HIC_CONSTANT_WATER_ACTIVITY’
model.root.input.model.unit_001.adsorption.is_kinetic = False # Kinetic binding
model.root.input.model.unit_001.adsorption.HICCWA_KA = [0, 3] # m^3 / (mol * s) (mobile phase)
model.root.input.model.unit_001.adsorption.HICCWA_KD = [0, 2000] # 1 / s (desorption)
model.root.input.model.unit_001.adsorption.HICCWA_NU = [0, 10]
model.root.input.model.unit_001.adsorption.HICCWA_QMAX = [0, 4]  # mol / m^3 
model.root.input.model.unit_001.adsorption.HICCWA_BETA0 = 0.28  # <- should be only one entry
model.root.input.model.unit_001.adsorption.HICCWA_BETA1 = 2.1e-4  # <- should be only one entry

HIC Binding WHS

model.root.input.model.unit_001.adsorption_model = ‘HIC_WATER_ON_HYDROPHOBIC_SURFACES’
model.root.input.model.unit_001.adsorption.is_kinetic = False # Kinetic binding
model.root.input.model.unit_001.adsorption.HICWHS_KA = [0, 0.3] # m^3 / (mol * s) (mobile phase)
model.root.input.model.unit_001.adsorption.HICWHS_KD = [0, 1.3] # 1 / s (desorption)
model.root.input.model.unit_001.adsorption.HICWHS_NU = [0, 10]
model.root.input.model.unit_001.adsorption.HICWHS_QMAX = [0, 4]  # mol / m^3 
model.root.input.model.unit_001.adsorption.HICWHS_BETA0 = 0.01  # <- should be only one entry
model.root.input.model.unit_001.adsorption.HICWHS_BETA1 = 1e-3  # <- should be only one entry

Let me know if that fixes the issues.

3 Likes