Colloidal isotherm hangs

CADET Troubleshooting
1

Expected behavior.

Simulation finishes with no binding. Intent is to fit Bpp and ln(Ke) at individual pH breakthrough and washout, then create pH dependence.

Actual behavior

Simulation hangs

How to produce bug (including a minimal reproducible example)

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# import cadet packages
from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import MultiComponentColloidal
from CADETProcess.processModel import Inlet, LumpedRateModelWithoutPores, Outlet
from CADETProcess.processModel import TubularReactor, Cstr
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process
from CADETProcess.reference import ReferenceIO
from CADETProcess.comparison import Comparator
from CADETProcess.optimization import OptimizationProblem
from CADETProcess.optimization import U_NSGA3
# --------------------------
# setup experiment
# --------------------------
mLmin_to_m3s = 60*1e6 # conversion
exp_flow = 8. # MV/min (mL/min)
Q = exp_flow/(mLmin_to_m3s) # convert to m^3/s
size_kDa = 150. # protein size kilodaltons
r_hyd_pro = 0.7429*size_kDa**0.3599 # protein hydrodynamic radius in m
c_salt = 1.0
# -------------------------------------
# load and wash and elute volumes and pH
# -------------------------------------
loadVol = 60.1 # mL
washVol = 35. # mL
eluteVol = 12. # mL
protein_c = 3.5 # g/L
pH_bind = 7.5
pH_elute = 3.8
pH_b_H = 10**(-pH_bind)
pH_e_H = 10**(-pH_elute)
# --------------------------------------
# calculating timepoints for experiment
# --------------------------------------
total_vol = loadVol + washVol + eluteVol # mL
exp_time = total_vol/exp_flow*60 # seconds
sample_vol = loadVol
load_time = loadVol/(exp_flow/60.)
wash_time = washVol/(exp_flow/60.)
elute_time = eluteVol/(exp_flow/60.)
# ---------------------------
# create component system
# ---------------------------
component_system = ComponentSystem()
component_system.add_component('salt')
component_system.add_component('pH')
component_system.add_component('protein')
# -------------------------
# binding model 
# -------------------------
binding_model = MultiComponentColloidal(component_system, name='MultiComponentColloidal')
# setup binding model
binding_model.is_kinetic = 1
binding_model.kinetic_rate_constant = [0,0,1]
binding_model.phase_ratio = 2.0e8
binding_model.coordination_number = 6
binding_model.protein_radius = [0,0,r_hyd_pro]
# kappa
binding_model.kappa_constant = 4e-9 # m debye param
binding_model.kappa_exponential = 0
binding_model.kappa_factor = 0
# logKeq params
binding_model.logkeq_ph_exponent = [0,0,0]
binding_model.logkeq_power_exponent = [0,0,0]
binding_model.logkeq_power_factor = [0,0,0]
binding_model.logkeq_exponent_factor = [0,0,0]
binding_model.logkeq_exponent_multiplier = [0,0,0]
# Bpp params
binding_model.bpp_ph_exponent = [0,0,0]
binding_model.bpp_power_exponent = [0,0,0]
binding_model.bpp_power_factor = [0,0,0]
binding_model.bpp_exponent_factor = [0,0,0]
binding_model.bpp_exponent_multiplier = [0,0,0]
# other params
binding_model.use_ph = False 
binding_model.linear_threshold = 1e-12 # should this be lower?
binding_model.phase_ratio = 2.0e8
binding_model.coordination_number = 6
binding_model.protein_radius = [0,0,r_hyd_pro]
binding_model.bound_states=[0,0,1]
print(binding_model.check_required_parameters())
# -------------------------
# setup BWE buffers
# -------------------------
# feed (mAb05)
feed = Inlet(component_system, name='feed')
feed.c = [c_salt,pH_b_H,protein_c/size_kDa]
# wash (washout buffer)
wash = Inlet(component_system, name='wash')
wash.c = [c_salt, pH_b_H,0]
# elute (elution buffer)
elute = Inlet(component_system, name='elute')
elute.c = [c_salt,pH_e_H,0]
# ----------------------
# device 
# ----------------------
device = LumpedRateModelWithoutPores(component_system, name='device')
device.length = 0.0015 # m 
device.diameter = 0.285 # m
device.total_porosity = 0.55
device.axial_dispersion = 1e-12 # low
device.binding_model = binding_model
device.discretization.ncol = 50
device.c = [c_salt,pH_b_H,0]
# --------------------------
# outlet
# --------------------------
outlet = Outlet(component_system, name='outlet')
outlet.solution_recorder.write_solution_bulk = True
#
# ------------------------------------------------
# now we create the Flow Sheet
# ------------------------------------------------
flow_sheet = FlowSheet(component_system)
flow_sheet.add_unit(feed)
flow_sheet.add_unit(wash)
flow_sheet.add_unit(elute)
flow_sheet.add_unit(device)
flow_sheet.add_unit(outlet)
# connect components
flow_sheet.add_connection(feed, device)
flow_sheet.add_connection(wash, device)
flow_sheet.add_connection(elute, device)
flow_sheet.add_connection(device,outlet)
# --------------------------------------------------------------
# now we create the process, i.e., turning feed on, swap to wash buffer, same flow rate
# --------------------------------------------------------------
process = Process(flow_sheet, 'breakthrough')
# bind
process.add_event('feed_on','flow_sheet.feed.flow_rate',Q,time=0)
process.add_event('feed_off','flow_sheet.feed.flow_rate',0,time=load_time)
# wash
process.add_event('wash_on','flow_sheet.wash.flow_rate',Q,time=load_time)
process.add_event('wash_off','flow_sheet.wash.flow_rate',0,time=load_time+wash_time)
# elute
#process.add_event('elute_on','flow_sheet.elute.flow_rate',Q,time=load_time + wash_time)
#process.add_event('elute_off','flow_sheet.elute.flow_rate',0,time=exp_time)
# cycle time
process.cycle_time = load_time + wash_time #exp_time # seconds
#process.n_cycles = 1
# -----------------------------
# simluate the process
# -----------------------------
if __name__ == '__main__':
    from CADETProcess.simulator import Cadet
    process_simulator = Cadet()
    process_simulator.time_resolution = 0.1
    simulation_results = process_simulator.simulate(process)

File produced by conda env export > environment.yml

environment.yml (7.56 KB)

3

Dear Rosario,

Thank you for reporting your issue!

I’ve identified two problems in your code:

  1. Your protein_radius (r_hyd_pro) calculates to approximately 4.5 m, which is unrealistically large and causes the simulation to break.
  2. Your bpp_exponent_factor is set to 0, which seems to cause issues with the solver. I’m not entirely sure why this happens yet, but for now, setting it to a small non-zero value (e.g., 1e-12) should resolve the problem. I’ll follow up here once we’ve identified the root cause.

With those two fixes implemented, your simulation has sucessfully run on my machine.

Since you’ve reported a few issues recently, I wanted to kindly remind you that it’s very helpful if you can provide a minimal reproducible example (MRE) when submitting bug reports. For example, in this case, elements like the events weren’t related to the problem. Stripping such unrelated parts out makes it quicker and easier for us to diagnose the issue and sometimes doing so can even help you spot the problem yourself. Thanks a lot!

All the best,
Hannah

3 Likes
4

Good morning Rosario,

As a follow-up to my previous post, I wanted to share a bit more detail on the issue you’re encountering.

The Multi-Component Colloidal binding model is defined as:

\begin{aligned} \frac{dq_i}{dt} &= k_{\text{kin},i} \left( c_{p,i} - c_{p,i}^\star \right), \quad i = 0, \dots, N_{\text{comp}} - 1 \\ c_{p,i}^\star &= q_i \exp \Biggl[ \frac{n(3 + \kappa R)}{4 q_{\text{tot}}} \sum_{j=0}^{N_{\text{bound}}} q_j \sqrt{b_{pp,i} b_{pp,j}} \frac{r_i + r_j}{2R} \exp \left( -\kappa \left( R - (r_i + r_j) \right) \right) - \ln(K_{e,i}) \Biggr] \end{aligned}

Here, the protein–protein interaction term b_{pp,i} is defined as:

b_{pp,i} = b_{a,i} \cdot c_{p,0}^{b_{b,i}} + b_{c,i} \cdot \exp\left( b_{d,i} \cdot c_{p,0} \right)

While setting b_{pp,i} = 0 would mathematically reduce the model to standard Langmuir binding kinetics, the solver is not equipped to handle zero vaules for b_{pp,i} . This is likely due to numerical issues such as division by zero or unexpectedly high or low derivatives in the Jacobian.

In your example, setting either bpp_power_factor b_{a,i} or bpp_exponent_factor b_{c,i} to zero leads to b_{pp,i} = 0 , which causes the solver to “hang”.

That said, since you want to fit b_{pp} , setting boundaries near zero instead of zero should keep the numerics stable.

Additionally, we have opened a GitHub issue to implement a fallback mechanism that catches the zero values and reduces the isotherm accordingly.

EDIT: Fixed in “Catches zero bpp value in Multi Component Colloidal #452”

Best regards,
Hannah

3 Likes
5

Thank you for assisting in my understanding behind the hanging simulations. quite interesting… In the future I’ll make sure to strip back to only essential code for a MRE.

Your explanation above and suggestions are certainly helpful.

Thanks again!

3 Likes
6

Good Morning Rosario,

You are very welcome, I am happy I was able to help!

All the best,
Hannah

1 Like