An example to use ACT isotherm for protein A affinity chromatography

Flynn · June 22, 2025, 6:10am

Introduction

The ACT isotherm is a modified Langmuir isotherm tailored for mechanistic modeling in protein A chromatography. It is known that both binding capacity q_{max} and equilibrium constant (affinity) K_{eq} changes as a function of pH. ACT isotherm uses modification factors derived based on proA and mAb interaction stoichiometries. Details can be found in the paper.

Example code

from cadet import Cadet
import numpy as np
import matplotlib.pyplot as plt

Cadet.cadet_path = r"C:\Users\zwend\CADET\cadet_act2025\bin\cadet-cli.exe"     ## IMPORTANT: use a version higher than 5.0.4 

## system settings
column_length = 0.025              ## m
column_volume = 5.025e-6           ## m^3
protein_MW = 150                   ## kDa
column_porosity = 0.31             ## 1
particle_porosity = 0.95           ## 1
Q = 1.25 /(6e7)                    ## volumetric flow rate, m^3/s

elution_pH_start = 5.5
elution_pH_end = 3.3

loading = 50.0                     ## mg/mL resin
c_feed = 4.0 / protein_MW          ## mol/m^3
RT = column_volume / Q             ## s

event_CV = [0.0, loading / (c_feed * protein_MW), 4.0, 10.0, ]           ## load, wash, elution in CV

Keq = 268                                                                ## ml/mg
Q_max = 85.6                                                             ## mg/ml column volume

total_porosity = column_porosity + (1.0-column_porosity)*particle_porosity

model = Cadet()

model.root.input.model.nunits = 3

model.root.input.model.unit_000.unit_type = 'INLET'
model.root.input.model.unit_000.ncomp = 2                            ## The first component is pH
model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'

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 = column_length                         # m              
model.root.input.model.unit_001.cross_section_area = column_volume/column_length   # m^2
model.root.input.model.unit_001.col_porosity = column_porosity                     # 1      
model.root.input.model.unit_001.par_porosity = particle_porosity                   # 1      
model.root.input.model.unit_001.par_radius = 0.0425e-3                             # m     
                                                                
## Transport
model.root.input.model.unit_001.col_dispersion = 1.36e-8            # m^2/s       
model.root.input.model.unit_001.film_diffusion = [1, 1.41e-5,]      # m/s
model.root.input.model.unit_001.par_diffusion = [1, 1.99e-11]       # m^2/s  
model.root.input.model.unit_001.par_surfdiffusion = [0.0, 0.0]      

model.root.input.model.unit_001.adsorption_model = 'AFFINITY_COMPLEX_TITRATION'

model.root.input.model.unit_001.adsorption.is_kinetic = 1
model.root.input.model.unit_001.adsorption.act_ka = [1.0, Keq*protein_MW*(1.0-total_porosity), ]          ##  m^3 solid phase / mol protein / s   ml/mg=m^3/kg    1 kda = 1 kg/mol
model.root.input.model.unit_001.adsorption.act_kd = [1.0, 1.0]                                            ## s^-1
model.root.input.model.unit_001.adsorption.act_qmax = [1e-10, Q_max/protein_MW/(1.0-total_porosity),]     ##  mol/m^3 solid phase mg/ml=kg/m^3 / 150 kg/mol = 1/150 mol/m^3

model.root.input.model.unit_001.adsorption.act_etaA = [0, 2.12,  ]
model.root.input.model.unit_001.adsorption.act_pkaA = [0, 3.82, ]
model.root.input.model.unit_001.adsorption.act_etaG = [0, 3.48, ]
model.root.input.model.unit_001.adsorption.act_pkaG = [0, 4.66, ]

model.root.input.model.unit_001.init_c = [7.2, 0.0, ]
model.root.input.model.unit_001.init_q = [0.0, 0.0]

### Grid cells
model.root.input.model.unit_001.discretization.ncol = 50
model.root.input.model.unit_001.discretization.npar = 12

### 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 = 2

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

model.root.input.solver.sections.nsec = 4
model.root.input.solver.sections.section_times = [0, event_CV[1] * RT, (event_CV[1] + event_CV[2]) * RT, (event_CV[1] + event_CV[2] + event_CV[3]) * RT, (event_CV[1] + event_CV[2] + event_CV[3] + 5.5) * RT]  ## s
model.root.input.solver.sections.section_continuity = [0,0,0,0]

## load
model.root.input.model.unit_000.sec_000.const_coeff = [7.2, c_feed,]                   # mol / m^3
model.root.input.model.unit_000.sec_000.lin_coeff =  [0.0, 0.0,]
model.root.input.model.unit_000.sec_000.quad_coeff = [0.0, 0.0,]
model.root.input.model.unit_000.sec_000.cube_coeff = [0.0, 0.0,]

## wash
model.root.input.model.unit_000.sec_001.const_coeff = [5.5, 0.0,]
model.root.input.model.unit_000.sec_001.lin_coeff =  [0.0, 0.0,]
model.root.input.model.unit_000.sec_001.quad_coeff = [0.0, 0.0,]
model.root.input.model.unit_000.sec_001.cube_coeff = [0.0, 0.0,]

## elution
model.root.input.model.unit_000.sec_002.const_coeff = [elution_pH_start, 0.0]
model.root.input.model.unit_000.sec_002.lin_coeff = [-(elution_pH_start-elution_pH_end)/(event_CV[3] * RT), 0.0, ]
model.root.input.model.unit_000.sec_002.quad_coeff = [0.0, 0.0,]
model.root.input.model.unit_000.sec_002.cube_coeff = [0.0, 0.0,]

## end
model.root.input.model.unit_000.sec_003.const_coeff = [elution_pH_end, 0.0]
model.root.input.model.unit_000.sec_003.lin_coeff = [0.0,0.0,]
model.root.input.model.unit_000.sec_003.quad_coeff = [0.0,0.0,]
model.root.input.model.unit_000.sec_003.cube_coeff = [0.0,0.0,]

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, Q,
    1, 2, -1, -1, Q]

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

model.root.input.solver.nthreads = 1

# Tolerances for the time integrator
model.root.input.solver.time_integrator.abstol = 1e-6
model.root.input.solver.time_integrator.algtol = 1e-10
model.root.input.solver.time_integrator.reltol = 1e-6
model.root.input.solver.time_integrator.init_step_size = 1e-6
model.root.input.solver.time_integrator.max_steps = 1000000

# Return data
model.root.input['return'].split_components_data = 1
model.root.input['return'].split_ports_data = 0
model.root.input['return'].unit_000.write_solution_bulk = 0
model.root.input['return'].unit_000.write_solution_inlet = 0
model.root.input['return'].unit_000.write_solution_outlet = 1

# Copy settings to the other unit operations
model.root.input['return'].unit_001 = model.root.input['return'].unit_000
model.root.input['return'].unit_002 = model.root.input['return'].unit_000

# Solution times
model.root.input.solver.user_solution_times = np.linspace(0, (event_CV[1] + event_CV[2] + event_CV[3] + 5.5) * RT, 201)

model.filename = 'model.h5'
model.save()
data = model.run()
model.load()   

time = model.root.output.solution.solution_times
c = model.root.output.solution.unit_001.solution_outlet_comp_001
pH_outlet = model.root.output.solution.unit_001.solution_outlet_comp_000
    
fig = plt.figure()
ax = fig.add_subplot()
ax.plot(time[1:]/60*Q*6e7, c[1:]*protein_MW, c="orange", label="mAb")
ax.set_xlabel(r'Volume/mL')
ax.set_ylabel(r'Concentration/(mg/mL)')

ax_ph = ax.twinx()
ax_ph.plot(time/60*Q*6e7, pH_outlet, label='pH')
ax_ph.set_ylabel('pH')
plt.show()

The above code should produce the following plot:

Feel free to change the code to step elution, different elution pHs etc.

Although the intended use of the ACT isotherm is to model pH effects, pH can be substituted for any type of salt as well. See the same paper Appendix for examples.

Please reach out to me with any questions regarding the isotherm usage.

naveen567 · July 1, 2025, 8:30pm

Hello Dr. Zhang,

Is it possible to use the ACT isotherm using CADET-Process? Also, I tried running the above example script using cadet-python, but it is not working. Kindly suggest.

Thank you!

Naveen J

Flynn · July 2, 2025, 4:41am

Hi Naveen,

Thank you for your interest. The above code requires a CADET version higher than 5.0.4, which is not yet released as of today. The team is working hard on a newer version, it will be available in near future. If you need it urgently, you can compile from source using the latest code on github. I can also send you a version I compiled alternatively.

To use the isotherm in CADET-Process, maybe @j.schmoelder can help?

j.schmoelder · July 2, 2025, 6:26am

Sure thing, it’s not a big issue to implement. Maybe you can join us tomorrow in our Office Hours @naveen567 and we can do it together?

naveen567 · July 2, 2025, 7:40am

Yeah Sure! We can do it tomorrow. Thank you!

naveen567 · July 2, 2025, 9:05am

Thank you, Dr. Zhang for the response. I will try compile the latest version using source code. Would you be able to provide the compiled version, if possible?

Thank you!

Naveen J

naveen567 · July 15, 2025, 6:33pm

Hi Johannes,

Quick update regarding the ACT isotherm through CADET-Process: the bound states are not defined.

class AffinityComplexTitration(BindingBaseClass):
    """
    Colloidal isotherm from Xu and Lenhoff 2009.

    Attributes
    ----------
    adsorption_rate: list of unsigned floats.
        Adsorption rate. Size depends on `n_comp`.
    desorption_rate: list of unsigned floats.
        Desorption rate. Size depends on `n_comp`.
    capacity: list of unsigned floats.
        Maximum adsorption capacities. Length depends on `n_comp`.
    eta_a: list of unsigned floats.
        Hill-type coefficients denoting the slope for the binding capacity changes as a
        function of pH changes. Length depends on `n_comp`.
    eta_g: list of unsigned floats.
        Hill-type coefficients denoting the slope for the equilibrium constant changes
        as a function of pH changes. Length depends on `n_comp`.
    pka_a: list of unsigned floats.
        Center point for the binding capacity changes as a function of pH changes.
        Length depends on `n_comp`.
    pka_g: list of unsigned floats.
        Center point for the equilibrium constant changes as a function of pH changes.
        Length depends on `n_comp`.
    """
    
    bound_states = SizedUnsignedIntegerList(
        size=("n_binding_sites", "n_comp"), default=1
    )

    adsorption_rate = SizedFloatList(size="n_comp")
    desorption_rate = SizedFloatList(size="n_comp")
    capacity = SizedUnsignedList(size="n_comp")
    eta_a = SizedUnsignedList(size="n_comp")
    eta_g = SizedUnsignedList(size="n_comp")
    pka_a = SizedUnsignedList(size="n_comp")
    pka_g = SizedUnsignedList(size="n_comp")

    _parameters = [
        "adsorption_rate",
        "desorption_rate",
        "capacity",
        "eta_a",
        "eta_g",
        "pka_a",
        "pka_g",
    ]

I edited the binding file locally, and it is working. So, I believe the remaining parts are fine.

Thank you!

Naveen J

j.schmoelder · July 16, 2025, 6:41am

Hi Naveen, thanks for reporting back.
If I understood correctly, you added those lines here:

    bound_states = SizedUnsignedIntegerList(
        size=("n_binding_sites", "n_comp"), default=1
    )

I don’t think we need this here since it should inherit this from the base class.
Could you please provide an example where this does not work?

Best

Jo

By the way, these technical aspects are best discussed over at the corresponding PR. Let us know if you have any questions.

naveen567 · July 16, 2025, 7:13pm

Hi all,

I am trying to use the ACT isotherm for protein A elution.

My expected behavior
Getting the elution profile during the elution step.

Actual behaviour

The elution keeps on going

Minimal reproducible example

import os
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from CADETProcess.processModel import ComponentSystem, Inlet, Outlet, Process, FlowSheet, GeneralRateModel, AffinityComplexTitration
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
from CADETProcess.simulator import Cadet

#install path
install_path = r"D:\Projects\CADET-Core\out\install\aRELEASE\bin\cadet-cli.exe"

mLmin_to_m3s = 60*1e6
exp_flow = 6
wash_flow = 9
Q = exp_flow/(mLmin_to_m3s)
Q2 = wash_flow/(mLmin_to_m3s)
size_kDa = 145
r_hyd_pro = 5e-9
protein_c = 5.3 #g/L
Q_max = 70   
#volume
loadvol = 540.38
wash1vol = 721.12 - 540.38
wash2vol = 780.76 - 721.12
wash3vol = 870.82 - 780.76
wash4vol = 1380.78667 - 870.82
elutevol = 1741.04 - 1380.78667
washVol = wash1vol+wash2vol+wash3vol+wash4vol
total_vol = loadvol + washVol + elutevol # mL
exp_time = 1741.04 # seconds
sample_vol = loadvol
load_time = 540.38
wash_time1 = 721.12
wash_time2 = 780.76
wash_time3 = 870.82
wash_time4 = 1380.79
elute_time = 1741.04

column_porosity =  0.6826             ## 1
particle_porosity = 0.168671867067481 
total_porosity = column_porosity + (1.0-column_porosity)*particle_porosity
#pH
pH_bind = 7
pH_elute = 3.5

#system component
component_system = ComponentSystem()
component_system.add_component('pH')
component_system.add_component('protein')

#binding model
binding_model = AffinityComplexTitration(component_system, name='ACT')
binding_model.bound_states = [0, 1]
binding_model.is_kinetic = 1
binding_model.adsorption_rate = [1, 1004.8281640948]          ##  m^3 solid phase / mol protein / s   ml/mg=m^3/kg    1 kda = 1 kg/mol
binding_model.desorption_rate = [1, 1.385]                                            ## s^-1
binding_model.capacity = [1e-10, 26.7119936577426]     ##  mol/m^3 solid phase mg/ml=kg/m^3 / 150 kg/mol = 1/150 mol/m^3

binding_model.eta_a = [0, 6.12]
binding_model.pka_a = [0, 3.82]
binding_model.eta_g = [0, 5.48]
binding_model.pka_g = [0, 6.66]
print(binding_model.check_required_parameters())

#buffers
# feed (mAb05)
feed = Inlet(component_system, name='feed')
feed.c = [pH_bind,protein_c/size_kDa]
# wash (washout buffer)
wash1 = Inlet(component_system, name='wash1')
wash1.c = [pH_bind,0]
wash2 = Inlet(component_system, name='wash2')
wash2.c = [pH_bind,0]
wash3 = Inlet(component_system, name='wash3')
wash3.c = [pH_bind,0]
wash4 = Inlet(component_system, name='wash4')
wash4.c = [pH_bind,0]
# elute (elution buffer)
elute = Inlet(component_system, name='elute')
elute.c = [pH_elute,0]

# Column
column = GeneralRateModel(component_system, name='column')
column.length = 0.12  # m 
column.diameter = 0.01  # m
column.particle_porosity = 0.6826
column.particle_radius = 50e-6
column.bed_porosity = 0.168671867067481
column.axial_dispersion = 2.74099487843477E-09
column.film_diffusion = [0, 0.0000048233133062222]
column.pore_diffusion = [0, 1.44342203138845E-11] 
column.binding_model = binding_model
column.discretization.ncol = 50
column.discretization.npar = 12
c = [7,protein_c/size_kDa]

# outlet
outlet = Outlet(component_system, name='outlet')
outlet.solution_recorder.write_solution_bulk = True
flow_sheet = FlowSheet(component_system)

flow_sheet.add_unit(feed, feed_inlet= True)
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, eluent_inlet=True)
flow_sheet.add_unit(column)
flow_sheet.add_unit(outlet, product_outlet=True)
flow_sheet.add_connection(feed, 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, 'Run1')
c_load = np.array([7, protein_c/size_kDa])
c_wash = np.array([7, 0])
t_gradient_start = wash_time4
gradient_duration = elute_time - wash_time4
gradient_slope = [(-3.5/(10)), 0]
c_gradient_poly = np.array(list(zip(c_wash, gradient_slope)))


#load
process.add_event('load', 'flow_sheet.feed.c', c_load)
process.add_event('load_on', 'flow_sheet.feed.flow_rate', Q, time=0)
process.add_event('load_off', 'flow_sheet.feed.flow_rate', 0, time=load_time)

#wash1
process.add_event('wash1', 'flow_sheet.wash1.c',  c_wash, load_time)
process.add_event('wash1_flow', 'flow_sheet.wash1.flow_rate', Q, time=load_time)
process.add_event('wash1_off', 'flow_sheet.wash1.flow_rate', 0, time=wash_time1)

#wash2
process.add_event('wash2', 'flow_sheet.wash2.c',  c_wash, wash_time1)
process.add_event('wash2_flow', 'flow_sheet.wash2.flow_rate', Q2, time=wash_time1)
process.add_event('wash2_off', 'flow_sheet.wash2.flow_rate', 0, time=wash_time2)

#wash3
process.add_event('wash3', 'flow_sheet.wash3.c',  c_wash, wash_time2)
process.add_event('wash3_on', 'flow_sheet.wash3.flow_rate', Q2, time=wash_time2)
process.add_event('wash3_off', 'flow_sheet.wash3.flow_rate', 0, time=wash_time3)

#wash4
process.add_event('wash4', 'flow_sheet.wash4.c',  c_wash, wash_time3)
process.add_event('wash4_on', 'flow_sheet.wash4.flow_rate', Q2, time=wash_time3)
process.add_event('wash4_off', 'flow_sheet.wash4.flow_rate', 0, time=wash_time4)

post_start = elute_time

#elute
process.add_event('elute_on', 'flow_sheet.elute.flow_rate', Q2, time=wash_time4)
#process.add_event('grad_start', 'flow_sheet.elute.c', c_gradient_poly, t_gradient_start)
process.add_event('grad_start', 'flow_sheet.elute.c', 3.5, t_gradient_start)
#process.add_event('grad_start_low', 'flow_sheet.elute.c', 3.5, time = elute_time-200 )
#process.add_event('grad_start', 'flow_sheet.elute.c', c_gradient_poly, t_gradient_start, polynomial_order=1)
process.add_event('elute_off', 'flow_sheet.elute.flow_rate', 0, time=elute_time)


# Cycle time
process.cycle_time = exp_time  # seconds

process_simulator = Cadet(install_path = install_path)
#simulator.evaluate_stationarity = True
process_simulator.check_cadet()

process_simulator.time_resolution = 0.1

simulation_results = process_simulator.simulate(process)
model_conc = simulation_results.solution.outlet.outlet.solution
model_times  = simulation_results.solution.outlet.outlet.time
from CADETProcess.plotting import SecondaryAxis
sec = SecondaryAxis()
sec.components = ['pH']
sec.y_label = '$c_{pH}$'

simulation_results.solution.column.outlet.plot(secondary_axis=sec)

CADET core used: v5.0.4
environment.yml (4.4 KB)

I tried gradient pH elution, but it is also not working. Can someone suggest what the error could be?

Thank you!

Naveeen J

Flynn · July 17, 2025, 4:45am

Hi Naveeen, unfortunately I am not familiar with cadet-process enough to do debugging, but if you can upload the associated .h5 file, I can take a look.

naveen567 · July 17, 2025, 8:53am

Hi Flynn,

Thank you for your response. I have attached the h5 file here.

ACT_elu.h5 (333.5 KB)

Naveen J