Reading this thread: Protein A affinity help and simulation error
It seems that the way to elution in affinity chromatography is to have the parameters of the binding model altered during elution. How can I model this in CADET-Process?
External functions are not yet supported in CADET-Process (see here).
1 Like
We have had good results with the pH-dependent colloidal isotherm that @soumi , @a.moser , and myself have worked with. I only have scripts for it in MATLAB form but they may be able to share python code. @Flynn also implemented a novel isotherm for proteinA that is available in some version of CADET. If you use either of these, you would not need to use external functions.
1 Like
That is correct, and @Ronan1234 I would be happy to discuss this with you or share code if you like. I have set up protein A models using the general rate model and the colloidal isotherm in both CADET (python) and CADET-process. I have added the colloidal isotherm and the [H+] concentration dependence for surface diffusivity to CADET-process for my own use, but it’s not part of the main version that’s available (yet).
Some notes on modeling protein A step elution:
Consistent with what @soumi showed in this paper, I have found that surface diffusivity is strongly dependent on [H+] concentration for protein A and is important for modeling protein A step elution. Additionally, it seems that it is also difficult to model protein A step elution without accounting for pH transients that can occur due to titration of the ligand and/or base matrix, though how important this is may be dependent on specific process conditions.
In short, having pH (or [H+]) dependent binding model parameters is necessary to model protein A elution, but it is only one piece of the puzzle.
2 Likes
Thank you for your response. Yes if you could share some code, that would be perfect!
Thanks for tagging me. Yes, I do have some MATLAB codes that I can share. Angela and @a4arjunvp have been taking the lead on python. 
Dear @soumi and @alter
Can you guys share the colloidal isotherm model code. If it is not in python script, you can share in MATLAB script. I am eagerly searching this code.
Any help is appreciated
Shyama
Hi Shyama,
@Daniel / @a4arjunvp should have the latest UDelaware versions.
Here is some relatively older MATLAB code (a function, and two different scripts, one to fit and one to predict). It should have the CADET inputs set up. I’ve invluded some sample data as well.
Thanks!
Soumi
Script_Predict_ColloidalLogEL.m (3.7 KB)
Script_SingleRun_ColloidalLogEL.m (3.9 KB)
formatted_1_MAbModel_MATLAB_NewCol_30MAR21_pH_0400 001.csv (483.8 KB)
ColloidalmodelEL_Rev_01.m (9.9 KB)
that I used.
Hi Soumi,
Thank you so much for sharing your code.
@Daniel and @a4arjunvp can you share latest versions of Colloidal isotherm model code.
Regards
Shyama
Here is my code, right now the code is adjusted for HIC but can very easily be adapted to ProA. I have included isotherm parameters from Soumi’s thesis. Also I am using a dev version of cadet process, not sure if the current release has the colloidal model yet.
fork: GitHub - angelamoser1/CADET-Process at add_GRM_surface_diffusion_dependence
If you need help on how to install it: A Guide to Reproducible Python Environments and CADET Installations - #8 by ronald.jaepel
post: " Installing a Custom Version of a Package From a Git Branch or Commit"
Best,
Daniel
Hi Daniel,
Thanks a lot for sharing the code. I have run the code in this dev branch for Protein, fork: GitHub - angelamoser1/CADET-Process at add_GRM_surface_diffusion_dependence. But getting an error: CADET Error: Simulation failed with b’IO ERROR: Field “COL_RADIUS” does not exist in group.
Please help to resolve the error.
Regards
Shyama
Hi Shyama,
That’s due to the change in the name of the variable from COL_RADIUS to COL_PROTEIN_RADIUS from CADET v4.4 to v5.0 that Ron mentioned here. The variable name in the CADET-Process branch is consistent with v5.0, COL_PROTEIN_RADIUS. You’d either need to upgrade CADET and be mindful of the kappa unit change, or go into the cadetAdapter file in CADET-Process and change the variable name back to COL_RADIUS to work with v4.4.
4 Likes
Hi Soumi,
I was going through your MATLAB code. For protein prediction script, you called a function PredictColloidalmodelEL_Rev_00(). Is it same function as ColloidalmodelEL_Rev_01() or different ?
If it is different, can you share PredictColloidalmodelEL_Rev_00() function. It is not attached what you shared.
Regards
Shyama
Hi Soumitra,
I used the load and wash BT data to fit the MCC parameters (Bpp, Lnke, and Keq). The fit result is as follows:
Further, I used the fitted parameters to simulate the full profile, but I am not getting any binding. I am simply getting a straight line.
The code I am using
import psutil
import collections
scpufreq = collections.namedtuple('scpufreq', ['current', 'min', 'max'])
def dummy_cpu_freq(percpu=False):
return scpufreq(current=2400.0, min=2400.0, max=2400.0)
psutil.cpu_freq = dummy_cpu_freq
import matplotlib.pyplot as plt
import pandas as pd
from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import MultiComponentColloidal
from CADETProcess.processModel import Inlet, Outlet, GeneralRateModel
from CADETProcess.processModel import FlowSheet, Process
from CADETProcess.simulator import Cadet
from CADETProcess.reference import ReferenceIO
from CADETProcess.comparison import Comparator
from CADETProcess.optimization import OptimizationProblem
from CADETProcess.optimization import U_NSGA3, Joblib
def set_binding_model(component_system, lnKe, Bpp, keq):
binding_model = MultiComponentColloidal(component_system, name='MultiComponentColloidal')
binding_model.is_kinetic = True
binding_model.phase_ratio = 1.46e8
binding_model.kappa_exponential = 0
binding_model.kappa_factor = 0
binding_model.kappa_constant = 4e-9
binding_model.coordination_number = 6
binding_model.logkeq_exponent_factor = [0, lnKe]
binding_model.logkeq_exponent_multiplier = [0, 0]
binding_model.logkeq_power_exponent = [0, 0]
binding_model.logkeq_power_factor = [0, 0]
binding_model.logkeq_ph_exponent = [0, 0]
binding_model.bpp_exponent_factor = [0, Bpp]
binding_model.bpp_exponent_multiplier = [0, 0]
binding_model.bpp_power_exponent = [0, 0]
binding_model.bpp_power_factor = [0, 0]
binding_model.bpp_ph_exponent = [0, 0]
binding_model.protein_radius = [0, 4.5e-9]
binding_model.kinetic_rate_constant = keq
binding_model.bound_states = [0, 1]
return binding_model
def set_column_model(component_system, binding_model):
column = GeneralRateModel(component_system, name='column')
column.length = 0.12 # m
column.diameter = 0.01 # m
column.particle_porosity = 0.168671867067481
column.particle_radius = 50e-6
column.bed_porosity = 0.6826
column.axial_dispersion = 4.38e-6
column.film_diffusion = [1, 9.10e-5]
column.pore_diffusion = [1, 1.25E-10]
column.binding_model = binding_model
column.discretization.ncol = 25
column.discretization.npar = 25
column.surface_diffusion = [0]
column.c = [1e3*10**(-7.0), 0]
column.cp = [1e3*10**(-7.0), 0]
column.q = [0]
return column
def make_process(lnKe, Bpp, keq):
component_system = ComponentSystem(['H+', 'protein'])
binding_model = set_binding_model(component_system, lnKe, Bpp, keq)
column = set_column_model(component_system, binding_model)
outlet = Outlet(component_system, name='outlet')
load = Inlet(component_system, name='load')
load.c = [1e3*10**(-7.0), 0.04]
wash1 = Inlet(component_system, name='wash1')
wash1.c = [1e3*10**(-7.0), 0]
wash2 = Inlet(component_system, name='wash2')
wash2.c = [1e3*10**(-7.0), 0]
wash3 = Inlet(component_system, name='wash3')
wash3.c = [1e3*10**(-7.0), 0]
wash4 = Inlet(component_system, name='wash4')
wash4.c = [1e3*10**(-7.0), 0]
elute = Inlet(component_system, name='elute')
elute.c = [1e3*10**(-3.5), 0]
flow_sheet = FlowSheet(component_system)
flow_sheet.add_unit(load)
flow_sheet.add_unit(wash1)
flow_sheet.add_unit(wash2)
flow_sheet.add_unit(wash3)
flow_sheet.add_unit(wash4)
flow_sheet.add_unit(elute)
flow_sheet.add_unit(column)
flow_sheet.add_unit(outlet)
flow_sheet.add_connection(load, column)
flow_sheet.add_connection(wash1, column)
flow_sheet.add_connection(wash2, column)
flow_sheet.add_connection(wash3, column)
flow_sheet.add_connection(wash4, column)
flow_sheet.add_connection(elute, column)
flow_sheet.add_connection(column, outlet)
process = Process(flow_sheet, name='proc')
load_time = 540.38
wash_time1 = 721.12
wash_time2 = 780.76
wash_time3 = 870.82
wash_time4 = 1380.79
elute_time = 1741.04
process.cycle_time = 1741.04
#events
process.add_event('load_on', 'flow_sheet.load.flow_rate', 1.00305555555556E-07, time=0)
process.add_event('load_off', 'flow_sheet.load.flow_rate', 0, time=load_time)
process.add_event('wash1_on', 'flow_sheet.wash1.flow_rate', 1.00305555555556E-07, time=load_time)
process.add_event('wash1_off', 'flow_sheet.wash1.flow_rate', 0, time=wash_time1)
process.add_event('wash2_on', 'flow_sheet.wash2.flow_rate', 1.50458333333333E-07, time=wash_time1)
process.add_event('wash2_off', 'flow_sheet.wash2.flow_rate', 0, time=wash_time2)
process.add_event('wash3_on', 'flow_sheet.wash3.flow_rate', 1.50458333333333E-07, time=wash_time2)
process.add_event('wash3_off', 'flow_sheet.wash3.flow_rate', 0, time=wash_time3)
process.add_event('wash4_on', 'flow_sheet.wash4.flow_rate', 1.50458333333333E-07, time=wash_time3)
process.add_event('wash4_off', 'flow_sheet.wash4.flow_rate', 0, time=wash_time4)
process.add_event('elute_on', 'flow_sheet.elute.flow_rate', 1.50458333333333E-07, time=wash_time4)
#process.add_event('grad_start', 'flow_sheet.elute.c', c_gradient_poly, t_gradient_start)
process.add_event('elute_off', 'flow_sheet.elute.flow_rate', 0, time=elute_time)
return process
def run_process(lnKe, Bpp, keq):
simulator = Cadet()
process = make_process(lnKe, Bpp, keq)
simulation_results = simulator.simulate(process)
final_time = simulation_results.solution.outlet.outlet.time
final_solution_array = simulation_results.solution.outlet.outlet.solution
conc_data = pd.read_csv('/Users/naveenj/Documents/Data/Data_elution.csv')
plt.plot(conc_data.iloc[:, 0], conc_data.iloc[:, 1], linestyle=':', color='k', label='reference')
plt.plot(final_time, final_solution_array[:, 1], color='b', label='fit')
plt.show()
return simulation_results
if __name__=='__main__':
run_process(17.83, 9.45, 10270609)
Kindly suggest.
Thanks
Naveen J
What you are calling k_{eq} is the kinetic rate constant k_{kin}, which you can fix at 1e8 and not fit this parameter. What are the values of lnKe and Bpp you are using? Also your values for mass transfer params are unrealistic. k_{film} and D_p should be on the order of 1e-4 and 1e-10 for pH, respectively, then 1e-5 and 1e-12 for mAb. I would recommend you take a look at my post on transport parameters estimation (below). Also, if you want to achieve a better fit to the front of the BT curve, it is recommended to fit D_p or obtain D_p from correlation and fit D_s (this is more physically consistent but also complicates things).
I also noticed the values for your porosities are very strange. Where are you getting these values? It actually matters a lot because it changes how the phase ratio is defined and significantly affects model parameter estimation.
3 Likes
Thanks @alters for sharing this.
Thank you for the suggestions @alters.
I am currently making modifications based on the suggestions. I modified my porosities; my new particle and bed porosities are 0.49 and 0.89, respectively.
1 Like
Hi Angela,
Can you please provide the code for cadet-python (parameter estimation for affinity chromatography colloidal model)?
Hi @alters,
I was just wondering if it is possible to model protein A elution using the Langmuir isotherm by tweaking the adsorption parameters. For example, putting low values of Keq and Qmax since during elution, there is weak binding. Also, putting the initial condition (q) of the column be the amount absorbed during loading instead of 0?