Hi, I am trying to replicate the HIC 2016 results both in GoSilico and Cadet. In cadet I was using the HIC Water on Hydrophobic surfaces isotherm and was doing this for 1 binding component so that the equations mentioned reduce to the original HIC 2016 model form. Upon setting the same configuration in both Cadet and Gosilico and using the same isotherm parameters with the right units I am getting a mismatch in the output chromatograms. I was using the HIC case study given in Gosilico for all my values. I am not sure if I am missing something in the equations or if the terms are defined differently in CADET. Has anyone else faced a similar issue? I have tried the same with the HIC 2008 Mollerup model and have found no discrepancies there.
Hey Shobhith,
that is to be expected, because the WHS isotherm is sensitive to the unit system used for calculations as reported in our recent publication here and GoSilico uses mol / L (as far as I know) and CADET uses SI units by default*, so mol / m³.
*: To be precise, CADET can be used with any unit system, but the user will need to ensure that all parameters and concentrations are given in matching units. To help with that, we label all parameters in SI units.
A possible work-around to get identical results is to use CADET in mol / L and load 0.0001 mol / L onto the column instead of 0.1 mol / m³. Adapt the salt gradient to mol / L, input all model parameters with the required transformations and you should get the same** results. If not, please let us know.
** small differences between model solvers can still remain, but the results should be similar.
Hi, Thanks for the reply. I have made sure to use the same units for both of them. In CADET I was using mol / m³ or mM and have converted the units from gosilico accordingly, but I am still getting varied results between the 2.
That is very interesting. Can you possibly share the CADET model setup code you used? I’d love to have a look.
Just a comment from the side line: it’s important which m³ you consider. Is it m³ of solid phase (which we assume in CADET) or particle volume (which is often assumed in other modeling frameworks).
case_study.ipynb (246.4 KB)
HIC2016.txt (278.5 KB)
I am attaching my Cadet simulation code as well the GoSilico XML file.
Ah, I see where a bit of the misunderstanding came from. It appears that GoSilico allows users to specify the units in which they enter the parameters. Under the hood (as far as I know), concentration parameters are then converted to mol / L. If I take your simulation and convert it to mol / L (by dividing all concentrations by 1000 and multiplying beta1 by 1000, as it is m^3/mol) I get something very similar to the GoSilico results:
I can’t reproduce your plots exactly because they refer to a ‘hic_10_pro3.xlsx’ file that I don’t have, but here’s the CADET side. However, the timing is still a bit off compared to your plot.
Taking into account what @j.schmoelder said, I tried to convert the parameters that use the volume of the stationary phase from lumped stationary phase to solid stationary phase, just to see what happens. Assuming GoSilico runs in \frac{mol}{V_{stationary}} and CADET in \frac{mol}{V_{solid}} with V_{stationary} = V_{solid} + V_{pores} and V_{pores} = V_{stationary}*\varepsilon_p we have a ratio of \frac{1}{1-\varepsilon_p}\frac{V_{stationary}}{V_{solid}}
We can now convert between the two with: \frac{mol}{V_{stationary}} = \frac{mol}{V_{stationary}}*\frac{1}{1-\varepsilon_p}\frac{V_{stationary}}{V_{solid}} = \frac{1}{1-\varepsilon_p}\frac{mol}{V_{solid}}
i.e. by dividing concentrations by (1-\varepsilon_p). This now gives profiles that have all elutions at the correct time, but with too wide elution profiles:
That GoSilico and CADET have different transport predictions is something we’ve seen before, but in the other direction, with CADET predicting too narrow of elutions.
Either way, as it is the weekend* and I’ve satiated my curiosity that the WHS isotherm can be approximately reproducible on both systems if the unit systems are considered, I’ll leave it at this.
*: (I’m not employed by the CADET Team anymore, so this is just me being curious in my free time)
Here’s the code for the plots in this post:
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
from CADETProcess.processModel import (
Inlet, Outlet, LumpedRateModelWithPores,
HICWaterOnHydrophobicSurfaces, ComponentSystem,
FlowSheet, Process
)
from CADETProcess.simulator import Cadet
warnings.filterwarnings("ignore")
concentration_factor = 1/1000
ref_vol_factor = 1/(1 - 0.7571)
def run_single_simulation(flow_sheet, Vg, Vl, tag):
simulator = Cadet()
process = create_process(flow_sheet, Vg, Vl, tag)
simulator.timeout = 6
result = simulator.simulate(process)
return tag, result, process
def create_process(flow_sheet, V_grad, V_load, tag):
process = Process(flow_sheet, f"{tag}")
wash_start = V_load / Q_m3_s
grad_start = wash_start + 5e-6 / Q_m3_s
grad_end = grad_start + V_grad / Q_m3_s
sim_state[f"grad_end_{tag}"] = grad_end
c_load = [starting_salt_conc, 0.136986*concentration_factor]
c_wash = [starting_salt_conc, 0.0]
c_elute = [75.2968*concentration_factor, 0.0]
slope = (np.array(c_elute) - np.array(c_wash)) / (grad_end - grad_start)
c_gradient_poly = list(zip(c_wash, slope))
process.cycle_time = grad_end + 5e-6 / Q_m3_s
process.add_event('load', 'flow_sheet.inlet.c', c_load)
process.add_event('wash', 'flow_sheet.inlet.c', c_wash, wash_start)
process.add_event('grad_start', 'flow_sheet.inlet.c', c_gradient_poly, grad_start)
process.add_event('elute_start', 'flow_sheet.inlet.c', c_elute, grad_end)
return process
component_system = ComponentSystem()
component_system.add_component('Salt')
component_system.add_component('Protein3')
binding = HICWaterOnHydrophobicSurfaces(component_system, name='HIC')
binding.is_kinetic = 1
equi = 3.89*ref_vol_factor
kinetic = 4.935
b0 = 0.0001048
b1 = 0.00450/concentration_factor
keff = 0.01e-3
n = 7.6
qmax = 48.7*concentration_factor*ref_vol_factor
Q_ml_min = 0.5
Q_m3_s = Q_ml_min / (60 * 1e6)
CV1 = 0.501437e-6
starting_salt_conc = 2083.45*concentration_factor
sim_state = {}
des_rate = 1 / kinetic
ads_rate = des_rate * equi
binding.adsorption_rate = [0.0, ads_rate]
binding.desorption_rate = [0.0, des_rate]
binding.hic_characteristic = [0.0, n]
binding.capacity = [0.0, qmax]
binding.beta_0 = b0
binding.beta_1 = b1
binding.bound_states = [0, 1]
inlet = Inlet(component_system, name='inlet')
inlet.flow_rate = Q_m3_s
column = LumpedRateModelWithPores(component_system, name='column')
column.binding_model = binding
column.length = 0.02
column.diameter = 0.00565
column.bed_porosity = 0.30
column.particle_radius = 4.5e-05
column.particle_porosity = 0.7571
column.axial_dispersion = 5.6e-8
column.film_diffusion = [0.01e-3, keff]
column.c = [starting_salt_conc, 0]
column.cp = [starting_salt_conc, 0]
column.q = [0]
column.discretization.ncol = 50
outlet = Outlet(component_system, name='outlet')
flow_sheet = FlowSheet(component_system)
flow_sheet.add_unit(inlet)
flow_sheet.add_unit(column)
flow_sheet.add_unit(outlet)
flow_sheet.add_connection(inlet, column)
flow_sheet.add_connection(column, outlet)
Vg, Vl, tag = 10 * CV1, 0.5e-6, 0
try:
tag, result, process = run_single_simulation(flow_sheet, Vg, Vl, tag)
sim_time = result.solution.outlet.outlet.time / 60 # Convert to minutes
sim_concentration = result.solution.outlet.outlet.total_concentration_components[:, 1]
sim_salt=result.solution.outlet.outlet.total_concentration_components[:, 0]
except Exception as e:
print("Simulation failed:", e)
# excel_file_path = 'hic_10_pro3.xlsx'
# excel_data = pd.read_excel(excel_file_path)
# time_from_excel = excel_data['Time [min]']
# concentration_from_excel = excel_data['Biomolecule Concentration [mM]']
fig, ax1 = plt.subplots(figsize=(10, 6))
ax1.plot(sim_time, sim_concentration, label="CADET", color='b')
# ax1.plot(time_from_excel, concentration_from_excel, label="GoSilico", color='r', linestyle='--')
ax1.set_xlabel('Time (minutes)')
ax1.set_ylabel('Protein (M)', color='b')
ax1.tick_params(axis='y', labelcolor='b')
plt.title('Simulation vs Excel Data')
ax1.legend(loc='upper left')
plt.show()
To summarize:
- The WHS isotherm is sensitive to the unit system used. i.e. mol / L will give different results than mol / m^3 → we need to run CADET in mol / L to reproduce GoSilico results (which always runs in mol / L) approximately.
Additionally, differences in reference volumes are probably a confounding factor. As the isotherm is unit-system-dependent, the difference in reference volumes is something we will never be able to compensate without re-writing our model code.As per Shobhiths answer, the change in 1. was enough.- Similar results are possible with both systems, but exact matches depend on too many factors that we don’t know about.
- Without access to a GoSilico license and lots of time to systematically test parameters (of which I have neither), I can’t say more, but I hope this helps.
Thank you very much for your reply. I have tried the code, and it works well. The timing error was due to a delay volume that was present in the Gosilico simulation, but otherwise the fix with just the unit system change was able to replicate the same results.
It is encouraging that we have some solid agreement between GoSilico and CADET since this has been an ongoing problem for multiple other examples.
One thing to note on your units is that I think solid phase in GoSilico and CADET is defined in terms of 1-\varepsilon_t not 1-\varepsilon_p since you need to consider the interstitial volume as well as intraparticle. This is what I have been doing in both simulators. In the HIC models I suppose this only affects q_{max,i} since we are not using ligand density.
P.S. I really enjoyed your recent paper!