An example to use ACT isotherm for protein A affinity chromatography

Hi Flynn,

When declaring values of Keq, qmax and kf in the script, do you use the values at binding pH or elution pH?

Nope, I am getting the same result. Perhaps the issue was that I only used the DBC load and wash data to calculate the binding parameters. Next, I plan to include an elution profile to enhance the binding parameters calculation.

According to my experience, the behavior you are seeing could be a result of combined isotherm and transport behaviors. First, ACT isotherm params can have an impact on the boarding of the peak, please refer to Figure 10 (b) 3 in my paper for an example. One thing you could try is to fit isotherm params closer or within to your loading density range to fine-tune it, see if it improves. Second, transport part is usually more complex, there might be contributions from both pore and surface diffusion. The boarding of the peak could just be because the diffusivity is decreased when q is too high, especially when they are released at the same time during elution, see this topic for some discussions. Another possibility could be surface diffusion, I also recommend playing around with it and see if it helps.

If you want to isolate these two, I recommend doing static binding experiments at different pHs to fix ACT isotherm params. And then in the model you can play around with transport params.

4 Likes

It should be binding pH, the elution is handled by the ACT params.

1 Like

For debugging, I would first try to use Langmuir isotherm and see if your code and most of the param settings are correct, make sure it is really binding, then move on to use ACT isotherm by adding ACT specific params for elution. After all, ACT isotherm is a modified Langmuir.

3 Likes

Hi Dr. Zhang,

As suggested, I have optimized the binding parameters using the NSGA-3, the binding parameters were optimized using load and wash data. I was able to get the following load+wash profile

However, when I was using the same parameters for full cycle, I am getting the below profile.

I guess the elution is not happening. Kindly suggest.

I have attached my h5 file for reference.

ACT_elu.h5 (219.9 KB)

Thank you!

Naveen J

The parameters for the pH is not set properly, for instance, I am seeing D_p and k_{film} are 0 for pH, which means it will not get into the pore at all. Please double check all the pH params (both transport and isotherm) against my example code.

2 Likes

Hi Flynn,

Thanks for your feedback. How do you find the transport and isotherm parameters for the pH component?

It is in my example code, pH is implemented as the first component as mentioned in my post as well.

It can be set to an arbitrarily small value, for example:

D_{p,pH} = 10^{-9}
k_{film,pH} = 10^{-4}

2 Likes

I tried to adjust the Langmuir parameters using the parameters given in the script I get this output

However, this is what I get from the actual values of Keq and Qmax obtained from SBC experiments

Can someone explain to me what is really an issue here?

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

Cadet.cadet_path = r"CADET-Core/install/bin/cadet-cli" 


## system settings
column_length = 0.025              ## m
column_volume = 1e-6          ## m^3
column_cross_section_area=column_volume/column_length #m2
protein_MW = 150                   ## kDa
column_porosity = 0.34           ## 1
particle_porosity = 0.94           ## 1
particle_radius= 150e-6       #m   actual give in COA =50um
# Keq = 3.22                    ## ml/mg
# Q_max = 37.35                 ## mg/ml resin
Col_dispersion =34.174e-6        ##m2/s 
film_diffusion = 1.96e-5       ##m2/s actual is 1.867e-6 affects positionn of peak if it is too low
total_porosity = column_porosity + (1.0-column_porosity)*particle_porosity
Q0 = 0.1 /(6e7)                    ## volumetric flow rate, m^3/s
Q1 = 2 /(6e7)                    ## volumetric flow rate, m^3/s


elution_pH_start = 3.8
elution_pH_end = 3.8

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

event_CV = [0.0, loading / (c_feed * protein_MW), 5.0, 5.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 = particle_radius                            # m     
                                                                
## Transport
model.root.input.model.unit_001.col_dispersion = Col_dispersion         # m^2/s       
model.root.input.model.unit_001.film_diffusion = [1, film_diffusion,]      # m/s
model.root.input.model.unit_001.par_diffusion = [1, 5.99e-8]       # 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, ]          ##  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,]     ##  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.4, 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.4, c_feed,]                   # mol / m^3

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

## elution
model.root.input.model.unit_000.sec_002.const_coeff = [elution_pH_start, 0.0]

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

model.root.input.model.connections.nswitches = 4
model.root.input.model.connections.switch_000.section = 0
model.root.input.model.connections.switch_000.connections = [
    0, 1, -1, -1, Q0,
    1, 2, -1, -1, Q0]

model.root.input.model.connections.switch_001.section = 1
model.root.input.model.connections.switch_001.connections = [
    0, 1, -1, -1, Q1,
    1, 2, -1, -1, Q1]

model.root.input.model.connections.switch_002.section = 2
model.root.input.model.connections.switch_002.connections = [
    0, 1, -1, -1, Q1,
    1, 2, -1, -1, Q1]

model.root.input.model.connections.switch_003.section = 3
model.root.input.model.connections.switch_003.connections = [
    0, 1, -1, -1, Q1,
    1, 2, -1, -1, Q1]


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, c[1:]*protein_MW, c="red", label="mAb")
ax.set_xlabel(r'Time (min)')
ax.set_ylabel(r'Concentration/(mg/mL)')

ax_ph = ax.twinx()
ax_ph.plot(time/60, pH_outlet, c="green", label='pH')
ax_ph.set_ylabel('pH')
plt.show()

I managed to reproduce the chromatogram as below after fixing some of the variables.




Thanks team for your support

I think the issue is that your Qmax is too small. The value you have from experiment is in terms of total column volume, it should be in terms of solid phase volume. Keq should also be adjusted to solid phase volume. You simply divide Qmax by (1-Et) and multiply Keq by (1-Et), also you need to convert to molar basis of course.

Thanks, @alters. The column that I am using has a very low Qmax, and from the SBC experiment, I got a very low value of Keq. May I know what is the acceptable range of these parameters?

What is the resin? For MS Sure for example, we would expect SBC around 30 mg/ml (on the low end) from batch adsorption and Keq of 100 ml/mg If your porosities are \varepsilon_e = 0.4 and \varepsilon_p = 0.9 this equates to \varepsilon_t = 0.94. You would then do the conversions

q_{max} = 30/(1-\varepsilon_t) = 500 \mathrm{\space mg / ml \space solid}

k_{eq} = 30\cdot(1-\varepsilon_t) = 6 \mathrm{\space ml \space solid / mg}

Then you convert these to molar basis (use mol/m^3, which is equivalent to mM) and use these values in CADET. The pH dependent adsorption params should be unitless but @Flynn would know better about this.

1 Like

Hi @alters,

I am using this resin column (ZetaSep Protein A Elevate).

Regarding ACT isotherm specific parameters, pKa_i is usually between 3 and 7 ish, while eta_i has a boarder range. I’ve seen values between [0, 15] for different mAbs on proA columns.

3 Likes

I don’t know much about this resin but it has a large particle size and seems sort of generic - I would use the parameters I mentioned as a starting point and look at some example parameters for MS SuRe.

1 Like

I also want to use the ACT model on the cadet-process. What should I do, please.

Hi,

The ACT isotherm is available for CADET-Core v5.1, which is available on conda forge. You can simply follow the installation instructions here. The interface for the ACT isotherm is documented here

@j.schmoelder can you say if it’s available with cadet-process as well?

Best regards,
Jan

1 Like