Trouble understanding the behaviour of SMA and getting the concentration peak for linear gradient elution

Chromatography Modeling
1

Hi, I was trying to use the SMA binding + Lumped Rate Model with Pores to model some chromatography system I needed to model.

I am having some trouble with these questions:

  1. If my salt concentration is 4 during both load and wash, why do I get that sharp peak of around 20 concentration for my salt?
  2. No matter what I do with absorption or desorption rate (even increasing and decreasing 100x, the chromatogram of protein 1 stays the exact same). Why is that?
  3. I have values for length, diameter, porosities, initial concentrations, load, wash, elute concentration etc. and the dispersion, diffusion, binding etc. parameters are unknown which I was trying to do a manual parameter estimation for (trying broad range of values and after getting a reasonable plot, fine-tune with optimisation). I am not getting good results at all, the entire protein is eluting before linear gradient elution even happens. The experimental data I have to work with has a narrow small peak during the end of load and start of wash and then a broad tall peak during elution.

P.S. The concentrations I had were in g/L not mM, I was intending to parameter estimate so that the same numbers can be used as values and the parameters will absorb any conversions that should have happened as I do not have the molecular weight for the protein.
PFA my code:

import numpy as np

from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import StericMassAction
from CADETProcess.processModel import Inlet, GeneralRateModel, Outlet, LumpedRateModelWithPores
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process

component_system = ComponentSystem()
component_system.add_component('Salt')
component_system.add_component('Protein 1')
component_system.add_component('Protein 2')

binding_model = StericMassAction(component_system, name='SMA')
binding_model.is_kinetic = True #False #True
binding_model.adsorption_rate = [0.0, 15000, 3] #[0.0, 30, 3]
binding_model.desorption_rate = [0.0, .015, 15] #[0.0, .5, 15]
binding_model.characteristic_charge = [0.0, 21.0, 7.0] #[0.0, 7.0, 7.0] #
binding_model.steric_factor = [0.0, 5, 50.0] #[0.0, 50.0, 50.0] #
binding_model.capacity = 37.64*3.565 #225.0 #

inlet = Inlet(component_system, name='inlet')
inlet.flow_rate = 1.63333e-8 #2.88e-8

# column = GeneralRateModel(component_system, name='column')
column = LumpedRateModelWithPores(component_system, name='column')
column.binding_model = binding_model

column.length = 0.2 #0.25
column.diameter = 0.005 #0.011286
column.bed_porosity = 0.39 #0.37
column.particle_radius = 2.5e-5 #4.5e-5
column.particle_porosity = 0.64 #0.33
column.axial_dispersion = 2.0e-7
column.film_diffusion = [2.0e-5, 2.0e-7, 2.0e-7]
# column.pore_diffusion = [7e-5, 1e-9, 1e-9]
# column.surface_diffusion = [0.0, 0.0, 0.0]

outlet = Outlet(component_system, name='outlet')

# Flow Sheet
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)

# Process
process = Process(flow_sheet, 'lwe')
process.cycle_time = 8160 #12000.0

## Create Events and Durations
wash_start = 2160 #4000.0
gradient_start = 3361 #8000.0
concentration_difference = np.array([37.5*3.565, 0.0]) - np.array([3.75*3.565, 0.0])
# concentration_difference = np.array([500, 0.0]) - np.array([70, 0.0])
gradient_duration = process.cycle_time - gradient_start
gradient_slope = concentration_difference/gradient_duration

_ = process.add_event('load', 'flow_sheet.inlet.c', [3.75*3.565, 6.86, 0.16]) #[180.0, 0.1, 0.1]) #
_ = process.add_event('wash', 'flow_sheet.inlet.c', [3.75*3.565, 0.0, 0.0], wash_start) #[70.0, 0.0, 0.0], wash_start) #
_ = process.add_event(
    'grad_start',
    'flow_sheet.inlet.c',
    [[3.75*3.565, gradient_slope[0]], [0, gradient_slope[1]], [0, gradient_slope[1]]],
    # [[70.0, gradient_slope[0]], [0, gradient_slope[1]], [0, gradient_slope[1]]],
    gradient_start
)

column.c = [3.75*3.565, 0, 0] #[180, 0, 0] #
column.q = [binding_model.capacity, 0, 0] #[3.75, 0, 0] #basetwo-eng-all 


Samardeep Sarna
  12:53 PM
import numpy as np

from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import StericMassAction
from CADETProcess.processModel import Inlet, GeneralRateModel, Outlet, LumpedRateModelWithPores
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process

component_system = ComponentSystem()
component_system.add_component('Salt')
component_system.add_component('Protein 1')
component_system.add_component('Protein 2')

binding_model = StericMassAction(component_system, name='SMA')
binding_model.is_kinetic = True 
binding_model.adsorption_rate = [0.0, 30, 3]
binding_model.desorption_rate = [0.0, .5, 15]
binding_model.characteristic_charge = [0.0, 21.0, 7.0] #[0.0, 7.0, 7.0] #
binding_model.steric_factor = [0.0, 5.0, 50.0] #[0.0, 50.0, 50.0] #
binding_model.capacity = 40 # g/L

inlet = Inlet(component_system, name='inlet')
inlet.flow_rate = 1.6e-8 

# column = GeneralRateModel(component_system, name='column')
column = LumpedRateModelWithPores(component_system, name='column')
column.binding_model = binding_model

column.length = 0.2 
column.diameter = 0.005 
column.bed_porosity = 0.4
column.particle_radius = 2.5e-5 
column.particle_porosity = 0.65
column.axial_dispersion = 2.0e-7
column.film_diffusion = [2.0e-5, 2.0e-7, 2.0e-7]

outlet = Outlet(component_system, name='outlet')

# Flow Sheet
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)

# Process
process = Process(flow_sheet, 'lwe')
process.cycle_time = 8200

## Create Events and Durations
wash_start = 2200
gradient_start = 3300
concentration_difference = np.array([40, 0.0]) - np.array([4, 0.0])
# concentration_difference = np.array([500, 0.0]) - np.array([70, 0.0])
gradient_duration = process.cycle_time - gradient_start
gradient_slope = concentration_difference/gradient_duration

_ = process.add_event('load', 'flow_sheet.inlet.c', [4, 7, 0.2]) #[180.0, 0.1, 0.1]) #
_ = process.add_event('wash', 'flow_sheet.inlet.c', [4, 0.0, 0.0], wash_start) #[70.0, 0.0, 0.0], wash_start) #
_ = process.add_event(
    'grad_start',
    'flow_sheet.inlet.c',
    [[4, gradient_slope[0]], [0, gradient_slope[1]], [0, gradient_slope[1]]],
    # [[70.0, gradient_slope[0]], [0, gradient_slope[1]], [0, gradient_slope[1]]],
    gradient_start
)

column.c = [4, 0, 0] 
column.q = [4, 0, 0] #[binding_model.capacity, 0, 0] #
2

Hi Ronan,

first: why does the code you shared contain two setups for a process? Which one should I work with? Please ensure the code you share is clear and concise in future.

Regarding your questions:

  1. If my salt concentration is 4 during both load and wash, why do I get that sharp peak of around 20 concentration for my salt?
  1. Because you displace the salt bound to the column during the loading step. Your protein fully binds to the column, displaces the salt ions, which then wash off of the column.
  1. No matter what I do with absorption or desorption rate (even increasing and decreasing 100x, the chromatogram of protein 1 stays the exact same). Why is that?
  1. Your code did not include the relevant section that does the plotting, so I can’t tell why the chromatogram looks identical for you. When I plot it with
from CADETProcess.simulator import Cadet

simulator = Cadet()
results = simulator.simulate(process)

from CADETProcess.plotting import SecondaryAxis

sec = SecondaryAxis()
sec.components = ['Salt']
sec.y_label = '$c_{salt}$'

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

I get different plots if I set the adsorption constants to 0, 100 or 100000. They are all, however, dominated by a breakthrough of both proteins during the loading phase because the capacity you set for your binding model does not have nearly enough room to hold the high concentrations of protein that you load for 36 minutes. To get realistic plots you need to lower the concentration (~0.1 works) or increase the capacity (I have not found a capacity high enough to handle the amount of protein loaded). Note: if you increase the capacity too much, simulations will get unstable unless you also introduce reference concentrations, so that would be recommended.

  1. I have values for length, diameter, porosities, initial concentrations, load, wash, elute concentration etc. and the dispersion, diffusion, binding etc. parameters are unknown which I was trying to do a manual parameter estimation for (trying broad range of values and after getting a reasonable plot, fine-tune with optimisation). I am not getting good results at all, the entire protein is eluting before linear gradient elution even happens. The experimental data I have to work with has a narrow small peak during the end of load and start of wash and then a broad tall peak during elution.
  1. Again, as above, that early elution is a break-through because you are overloading your column. Another question: Are you trying to estimate transport and binding parameters in one step off of gradient elution data? That is unlikely to give good results and will take a substantial amount of time because you have a high dimensional problem with colinearities between parameters. Usually it is better to first estimate transport parameters on non-binding tracer experiments (e.g. your protein under high-salt non-binding conditions) to isolate the transport effects and then move to binding parameters once you have the transport parameters.

P.S. The concentrations I had were in g/L not mM, I was intending to parameter estimate so that the same numbers can be used as values and the parameters will absorb any conversions that should have happened as I do not have the molecular weight for the protein.
PFA my code:

The SMA model does not work with mixed units (mM for salt and g/L for protein) and I don’t understand how it works if it’s entirely in g/L either. It seems to me to break the electro-neutrality conditional calculations, but I’m not sure. Maybe it can work. I would still highly recommend getting an approximate estimate of the molecular weight.

2 Likes