Background
As I discussed in the Recovery conference, the design philosophy of ACT isotherm goes far beyond affinity chromatography. In theory, it can be applied to any modes of chromatography as detailed in the paper Appendix section (generalized ACT or gACT). Langmuir isotherms can be generally applied to any modes of chromatography including CEX, gACT modifies Langmuir isotherm and provides general handles like pH, H+, Na+, Ca2+, etc. to impact the binding capacity and equilibrium constants.
Recently the implementation of ACT isotherm has been improved to accommodate gACT. Below are some examples to use gACT to simulate a linear salt gradient in a cation exchange chromatography step.
Split peaks
Split peaks may arise from many different sources. The ACT isotherm can mimic some of these effects because 1. if the apparent binding capacity and apparent equilibrium constant change too differently, it tends to give shoulders or split peaks; 2. ACT assumes competitive binding to the same binding sites for different components, which could react to pH/salt concentration dramatically differently during elution. The more strongly binding species can delay, sharpen, or partially displace the weaker species, producing apparent split peak or asymmetric elution profiles.
With the above said, those are perspectives from the isotherm. Whether or not it makes physical sense or what is causing these behaviors physically depends on the users’ interpretation.
Single component system w/ and w/o split peaks
A simple single component system:
Change the isotherm parameter we could get split peaks even for a single component system:
Multi component system w/ split peaks
For a multicomponent system, one can also get split peaks. For instance:
Monomer and LMW12 are configured to have only one peak. Zoomed-in versions for HMW13 and HMW2 which have split peaks:
Depending on the parameters, there could be more than 2 split peaks. I will share a good example when I come across one.
Code
The CADET-Core code will be released along with version 6 in the future. For urgent use, please compile the latest version of CADET-Core on github or DM me. Below is the CADET-Python code to reproduce the above figures.
Single component system:
Cadet.cadet_path = ## use your cadet_cli.exe path
def run_cadet(x, start_salt_conc, end_salt_conc, grad_start_time, grad_end_time, loaded_mass, filename):
'''
@input: x contains all input params for monomer only
x = [etaA, cmidA, etaG, cmidG, log(Dp)]
loaded_mass in mg (scalar, monomer only)
'''
etaA = x[0]
cmidA = x[1]
etaG = x[2]
cmidG = x[3]
Dp = np.exp(x[4])
Q_max = 104.91
Keq = 1189.96
n_comp = 1 # monomer only
protein_MW = 150e3 # Da = g/mol = mg/mmol
load_conc = 5.99 # mg/ml
column_porosity = 0.40
particle_porosity = 0.65
RT = 3 * 60 # s
total_porosity = column_porosity + (1.0 - column_porosity) * particle_porosity
model = Cadet()
## Number of units
model.root.input.model.nunits = 3
## Inlet
model.root.input.model.unit_000.unit_type = 'INLET'
model.root.input.model.unit_000.ncomp = 1 + n_comp # salt + monomer
model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
## Column model (GRM)
model.root.input.model.unit_001.unit_type = 'COLUMN_MODEL_1D'
model.root.input.model.unit_001.ncomp = 1 + n_comp
model.root.input.model.unit_001.npartype = 1
## Geometry
model.root.input.model.unit_001.col_length = 0.1
model.root.input.model.unit_001.cross_section_area = 0.25 * np.pi * 0.005**2
CV = model.root.input.model.unit_001.col_length * model.root.input.model.unit_001.cross_section_area
Q = CV / RT
loading = loaded_mass / (CV * 1e6) # mg/ml resin
model.root.input.model.unit_001.col_porosity = column_porosity
model.root.input.model.unit_001.particle_type_000.par_porosity = particle_porosity
model.root.input.model.unit_001.particle_type_000.par_radius = 50.0 / 2.0 * 1e-6
## Transport
model.root.input.model.unit_001.col_dispersion = 6e-8
model.root.input.model.unit_001.particle_type_000.has_film_diffusion = True
model.root.input.model.unit_001.particle_type_000.film_diffusion = np.array([1e-05, 1e-05], dtype=np.float64)
model.root.input.model.unit_001.particle_type_000.has_pore_diffusion = True
model.root.input.model.unit_001.particle_type_000.pore_diffusion = np.array([7e-10, Dp], dtype=np.float64)
model.root.input.model.unit_001.particle_type_000.has_surface_diffusion = False
## Isotherm
model.root.input.model.unit_001.particle_type_000.adsorption_model = 'AFFINITY_COMPLEX_TITRATION'
model.root.input.model.unit_001.particle_type_000.nbound = np.array([0, 1])
model.root.input.model.unit_001.particle_type_000.adsorption.act_use_ion_conc = True
model.root.input.model.unit_001.particle_type_000.adsorption.is_kinetic = 1
model.root.input.model.unit_001.particle_type_000.adsorption.act_ka = [
1.0,
Keq * protein_MW * (1.0 - total_porosity),
]
model.root.input.model.unit_001.particle_type_000.adsorption.act_kd = [1.0, 1.0]
model.root.input.model.unit_001.particle_type_000.adsorption.act_qmax = [
1e-10,
Q_max / protein_MW / (1.0 - total_porosity),
]
model.root.input.model.unit_001.particle_type_000.adsorption.act_etaa = np.array([0.0, etaA], dtype=np.float64)
model.root.input.model.unit_001.particle_type_000.adsorption.act_etag = np.array([0.0, etaG], dtype=np.float64)
model.root.input.model.unit_001.particle_type_000.adsorption.act_cmid_a = np.array([0.0, cmidA], dtype=np.float64)
model.root.input.model.unit_001.particle_type_000.adsorption.act_cmid_g = np.array([0.0, cmidG], dtype=np.float64)
## Initial conditions
model.root.input.model.unit_001.init_c = [start_salt_conc, 0.0]
model.root.input.model.unit_001.init_cp = [start_salt_conc, 0.0]
## FVM discretization
model.root.input.model.unit_001.discretization.spatial_method = "FV"
model.root.input.model.unit_001.discretization.reconstruction = np.bytes_(b'WENO')
model.root.input.model.unit_001.discretization.ncol = np.int64(40)
model.root.input.model.unit_001.discretization.weno.weno_order = np.int64(3)
model.root.input.model.unit_001.particle_type_000.discretization.spatial_method = "FV"
model.root.input.model.unit_001.particle_type_000.discretization.ncells = np.int64(20)
model.root.input.model.unit_001.particle_type_000.discretization.par_disc_type = np.bytes_(b'EQUIDISTANT_PAR')
model.root.input.model.unit_001.discretization.use_analytic_jacobian = 1
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.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
## Time integrator tolerances
model.root.input.solver.time_integrator.abstol = 1e-08
model.root.input.solver.time_integrator.algtol = 1e-12
model.root.input.solver.time_integrator.init_step_size = 1e-12
model.root.input.solver.time_integrator.max_convtest_fail = 50
model.root.input.solver.time_integrator.max_errtest_fail = 7
model.root.input.solver.time_integrator.max_steps = 1000000
model.root.input.solver.time_integrator.reltol = 1e-08
## Outlet
model.root.input.model.unit_002.unit_type = 'OUTLET'
model.root.input.model.unit_002.ncomp = 1 + n_comp
## Section times: load, wash, elution, hold
lag_time = 1300
model.root.input.solver.sections.nsec = 4
model.root.input.solver.sections.section_times = [
0,
RT * loading / load_conc,
grad_start_time,
grad_end_time,
grad_end_time + lag_time,
]
model.root.input.solver.sections.section_continuity = [0, 0, 0, 0]
## Load phase
model.root.input.model.unit_000.sec_000.const_coeff = [start_salt_conc, load_conc / protein_MW]
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 phase
model.root.input.model.unit_000.sec_001.const_coeff = [start_salt_conc, 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 phase
model.root.input.model.unit_000.sec_002.const_coeff = [start_salt_conc, 0.0]
model.root.input.model.unit_000.sec_002.lin_coeff = [
(end_salt_conc - start_salt_conc) / (grad_end_time - grad_start_time), 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]
## Hold phase
model.root.input.model.unit_000.sec_003.const_coeff = [end_salt_conc, 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]
## Return / output settings
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
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, grad_end_time + lag_time, 500)
## Connections
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.filename = filename
model.save()
data = model.run()
model.load()
time = model.root.output.solution.solution_times
volume = time * Q * 1e6 # mL
salt = model.root.output.solution.unit_001.solution_outlet_comp_000
monomer = model.root.output.solution.unit_001.solution_outlet_comp_001 * protein_MW
return volume, salt, monomer
## sinlge component, simulation 1, single peak
x = np.asarray([
2.28, # etaA
270.62, # cmidA
10.89, # etaG,
310.76, # cmidG,
-26.3, # log(Dp)
])
loaded_mass = 53.12209550021964 # mg, monomer only
volume, salt, monomer = run_cadet(x, 54, 1038, 2373, 9647, loaded_mass, 'model.h5')
## plot
fig, axes = plt.subplots(1, 1, figsize=(12, 6), )
axes.plot(volume, monomer, color='red', label='Monomer', alpha=0.6)
axes.set_ylabel('Concentration (mg/mL)')
axes.set_title('Single component, linear salt gradient using gACT isotherm')
axes.legend(loc = 'upper left')
axes.grid(True, alpha=0.3)
axes.set_xlabel('Volume (mL)')
axes2 = axes.twinx()
axes2.plot(volume, salt, color='blue', label='Na⁺ (salt)', alpha=0.6)
axes2.set_ylabel('Concentration (mM)')
axes2.legend(loc = 'center left')
plt.show()
## sinlge component, simulation 2, split peaks
x = np.asarray([
2.28, # etaA
220.62, # cmidA
13.89, # etaG,
400.76, # cmidG, 310.76
-26.3, # log(Dp)
])
loaded_mass = 53.12209550021964 # mg, monomer only
volume, salt, monomer = run_cadet(x, 54, 1038, 2373, 9647, loaded_mass, 'model.h5')
# plot
fig, axes = plt.subplots(1, 1, figsize=(12, 6), )
axes.plot(volume, monomer, color='red', label='Monomer', alpha=0.6)
axes.set_ylabel('Concentration (mg/mL)')
axes.set_title('Single component with split peaks, linear salt gradient using gACT isotherm')
axes.legend(loc = 'upper left')
axes.grid(True, alpha=0.3)
axes.set_xlabel('Volume (mL)')
axes2 = axes.twinx()
axes2.plot(volume, salt, color='blue', label='Na⁺ (salt)', alpha=0.6)
axes2.set_ylabel('Concentration (mM)')
axes2.legend(loc = 'center left')
axes2.set_xlabel('Volume (mL)')
plt.show()
Multicomponent system:
def run_cadet(x, start_salt_conc, end_salt_conc, grad_start_time, grad_end_time, loaded_mass, filename):
'''
@input: x contains all input params: the order should be qmax, keq, etaA, pkaA, etaG, pkaG
loaded_mass in mg, mass of all components
'''
## comfigure for a multicomponent system
## the sequence is HMW31, HMW2, mono, LMW12
etaA = [x[0], x[1], x[2], x[3], ]
cmidA = [x[4], x[5], x[6], x[7], ]
etaG = [x[8], x[9], x[10], x[11], ]
cmidG = [x[12], x[13], x[14], x[15], ]
Dps = np.exp(np.asarray([x[16], x[17], x[18], x[19], ]))
Q_max = [188.96, 197.45, 104.91, 112.26, ]
Keq = [1253.32, 1057.12, 1189.96, 2177.84, ]
n_comp = len(etaA)
protein_MW = np.asarray([450e3, 300e3, 150e3, 75e3]) ## Da = g/mol = mg/mmol
load_conc = 5.99 ## mg/ml
column_porosity = 0.40 ## 1
particle_porosity = 0.65 ## 1
RT = 3*60 ## s
total_porosity = column_porosity + (1.0-column_porosity)*particle_porosity
load_ratio = [loaded_mass[i] / np.sum(loaded_mass) for i in range (0, n_comp)]
model = Cadet()
## number of units
model.root.input.model.nunits = 3
## Inlet
model.root.input.model.unit_000.unit_type = 'INLET'
model.root.input.model.unit_000.ncomp = 1+n_comp # salt + proteins
model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
## column model
## GRM
model.root.input.model.unit_001.unit_type = 'COLUMN_MODEL_1D'
model.root.input.model.unit_001.ncomp = 1+n_comp
## Geometry
model.root.input.model.unit_001.col_length = 0.1 # m
model.root.input.model.unit_001.cross_section_area = 0.25*np.pi*0.005**2 # m
CV = model.root.input.model.unit_001.col_length * model.root.input.model.unit_001.cross_section_area ## m^3
Q = CV / RT # m^3 / s
loading = np.sum(loaded_mass) / (CV*1e6) # mg/ml resin
model.root.input.model.unit_001.col_porosity = column_porosity # -
model.root.input.model.unit_001.particle_type_000.par_porosity = particle_porosity # -
model.root.input.model.unit_001.particle_type_000.par_radius = 50.0/2.0*1e-6 # m
model.root.input.model.unit_001.npartype=1
## Transport
model.root.input.model.unit_001.col_dispersion = 6e-8 # m^2 / s (interstitial volume)
model.root.input.model.unit_001.particle_type_000.has_film_diffusion = True
model.root.input.model.unit_001.particle_type_000.film_diffusion = np.array([1e-05, 1e-05, 1e-05, 1e-05, 1e-05, ], dtype=np.float64) # m / s
model.root.input.model.unit_001.particle_type_000.has_pore_diffusion = True
model.root.input.model.unit_001.particle_type_000.pore_diffusion = np.array([7e-10, *Dps], dtype=np.float64)
model.root.input.model.unit_001.particle_type_000.has_surface_diffusion = False
## Isotherm model
model.root.input.model.unit_001.particle_type_000.adsorption_model = 'AFFINITY_COMPLEX_TITRATION'
model.root.input.model.unit_001.particle_type_000.nbound = np.array([0, 1, 1, 1, 1])
model.root.input.model.unit_001.particle_type_000.adsorption.act_use_ion_conc = True
model.root.input.model.unit_001.particle_type_000.adsorption.is_kinetic = 1
model.root.input.model.unit_001.particle_type_000.adsorption.act_ka = [1.0,
Keq[0]*protein_MW[0]*(1.0-total_porosity),
Keq[1]*protein_MW[1]*(1.0-total_porosity),
Keq[2]*protein_MW[2]*(1.0-total_porosity),
Keq[3]*protein_MW[3]*(1.0-total_porosity), ] ## ml liq/mg * mg/mmol * mL resin/mL liq = ml resin/mmol (mM-1)
model.root.input.model.unit_001.particle_type_000.adsorption.act_kd = [1.0, 1.0, 1.0, 1.0, 1.0, ]
model.root.input.model.unit_001.particle_type_000.adsorption.act_qmax = [1e-10,
Q_max[0]/protein_MW[0]/(1.0-total_porosity),
Q_max[1]/protein_MW[1]/(1.0-total_porosity),
Q_max[2]/protein_MW[2]/(1.0-total_porosity),
Q_max[3]/protein_MW[3]/(1.0-total_porosity),] ## mmol/ml resin (mM)
model.root.input.model.unit_001.particle_type_000.adsorption.act_etaa = np.array([0.0, *etaA], dtype=np.float64)
model.root.input.model.unit_001.particle_type_000.adsorption.act_etag = np.array([0.0, *etaG], dtype=np.float64)
model.root.input.model.unit_001.particle_type_000.adsorption.act_cmid_a = np.array([0.0, *cmidA], dtype=np.float64) ## in mM
model.root.input.model.unit_001.particle_type_000.adsorption.act_cmid_g = np.array([0.0, *cmidG], dtype=np.float64) ## in mM
model.root.input.model.unit_001.init_c = [start_salt_conc, 0, 0, 0, 0, ]
model.root.input.model.unit_001.init_cp = [start_salt_conc, 0, 0, 0, 0, ]
## FVM discretization
### Grid cells
model.root.input.model.unit_001.discretization.spatial_method = "FV"
model.root.input.model.unit_001.discretization.reconstruction = np.bytes_(b'WENO')
model.root.input.model.unit_001.discretization.ncol = np.int64(40)
model.root.input.model.unit_001.discretization.weno.weno_order = np.int64(3)
model.root.input.model.unit_001.particle_type_000.discretization.spatial_method = "FV"
model.root.input.model.unit_001.particle_type_000.discretization.ncells = np.int64(20)
model.root.input.model.unit_001.particle_type_000.discretization.par_disc_type = np.bytes_(b'EQUIDISTANT_PAR')
### Other options
model.root.input.model.unit_001.discretization.use_analytic_jacobian = 1
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.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
# Number of cores for parallel simulation
model.root.input.solver.nthreads = 1
# Tolerances for the time integrator
model.root.input.solver.time_integrator.abstol = 1e-08
model.root.input.solver.time_integrator.algtol = 1e-12
model.root.input.solver.time_integrator.init_step_size = 1e-12
model.root.input.solver.time_integrator.max_convtest_fail = 50
model.root.input.solver.time_integrator.max_errtest_fail = 7
model.root.input.solver.time_integrator.max_steps = 1000000
model.root.input.solver.time_integrator.reltol = 1e-08
## outlet
model.root.input.model.unit_002.unit_type = 'OUTLET'
model.root.input.model.unit_002.ncomp = 1+n_comp
## time
model.root.input.solver.sections.nsec = 4
# load, wash, elution, hold
lag_time = 1300
model.root.input.solver.sections.section_times = [0, RT*loading/load_conc, grad_start_time, grad_end_time, grad_end_time+lag_time, ] # s
model.root.input.solver.sections.section_continuity = [0, 0, 0, 0,]
## load phase
model.root.input.model.unit_000.sec_000.const_coeff = [start_salt_conc,
load_conc * load_ratio[0] /protein_MW[0],
load_conc * load_ratio[1] /protein_MW[1],
load_conc * load_ratio[2] /protein_MW[2],
load_conc * load_ratio[3] /protein_MW[3], ] # pksalt and M: mg/mL * mmol/mg = mmol/ml
model.root.input.model.unit_000.sec_000.lin_coeff = 5*[0.0,]
model.root.input.model.unit_000.sec_000.quad_coeff = 5*[0.0,]
model.root.input.model.unit_000.sec_000.cube_coeff = 5*[0.0,]
## wash phase
model.root.input.model.unit_000.sec_001.const_coeff = [start_salt_conc, 0.0, 0.0, 0.0, 0.0, ] # mol / m^3
model.root.input.model.unit_000.sec_001.lin_coeff = 5*[0.0,]
model.root.input.model.unit_000.sec_001.quad_coeff = 5*[0.0,]
model.root.input.model.unit_000.sec_001.cube_coeff = 5*[0.0,]
## elution phase
model.root.input.model.unit_000.sec_002.const_coeff = [start_salt_conc, 0.0, 0.0, 0.0, 0.0, ] # mol / m^3
model.root.input.model.unit_000.sec_002.lin_coeff = [(end_salt_conc - start_salt_conc) / (grad_end_time - grad_start_time), 0.0, 0.0, 0.0, 0.0, ]
model.root.input.model.unit_000.sec_002.quad_coeff = 5*[0.0,]
model.root.input.model.unit_000.sec_002.cube_coeff = 5*[0.0,]
## hold phase
model.root.input.model.unit_000.sec_003.const_coeff = [end_salt_conc, 0, 0, 0, 0, ] # mol / m^3
model.root.input.model.unit_000.sec_003.lin_coeff = 5*[0.0,]
model.root.input.model.unit_000.sec_003.quad_coeff = 5*[0.0,]
model.root.input.model.unit_000.sec_003.cube_coeff = 5*[0.0,]
# 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 # off, spatial concentration distributions
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, grad_end_time+lag_time, 500)
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, # [unit_000, unit_001, all components, all components, Q/ m^3*s^-1
1, 2, -1, -1, Q, ] # [unit_001, unit_002, all components, all components, Q/ m^3*s^-1
model.filename = filename
model.save()
data = model.run()
#print(data)
model.load()
time = model.root.output.solution.solution_times
volume = time * Q * 1e6 # mL
salt = model.root.output.solution.unit_001.solution_outlet_comp_000
HMW13 = model.root.output.solution.unit_001.solution_outlet_comp_001 * protein_MW[0]
HMW2 = model.root.output.solution.unit_001.solution_outlet_comp_002 * protein_MW[1]
monomer = model.root.output.solution.unit_001.solution_outlet_comp_003 * protein_MW[2]
LMW12 = model.root.output.solution.unit_001.solution_outlet_comp_004 * protein_MW[3]
return volume, salt, HMW13, HMW2, monomer, LMW12
## HMW13, HMW2, monomer, LMW12
x = np.asarray([ 7.07, 8.97, 12.28, 5.40, ## etaA
460.53, 400.50, 260.62, 230.19, ## cmidA
15.27, 20.64, 31.82, 16.64, ## etaG
186.84, 240.46, 255.74, 170.96, ## cmidG
-27.7, -27.3, -26.9, -26.6, ]) ## log(Dp)
loaded_mass = [np.float64(2.9554704554706563),
np.float64(2.757192885678953),
np.float64(53.12209550021964),
np.float64(10.436484908630728)]
volume, salt, HMW13, HMW2, monomer, LMW12 = run_cadet(x, 54, 1038, 2373, 9647, loaded_mass, 'model.h5')
fig, axes = plt.subplots(1, 1, figsize=(12, 6))
# Primary y-axis: protein species
axes.plot(volume, monomer, color='red', label='Monomer', alpha=0.6)
axes.plot(volume, HMW13, color='green', label='HMW13', alpha=0.6)
axes.plot(volume, HMW2, color='orange', label='HMW2', alpha=0.6)
axes.plot(volume, LMW12, color='purple', label='LMW12', alpha=0.6)
axes.set_ylabel('Concentration (mg/mL)')
axes.set_title('Multi-component, linear salt gradient using gACT isotherm')
axes.legend(loc='upper left')
axes.grid(True, alpha=0.3)
axes.set_xlabel('Volume (mL)')
# Secondary y-axis: salt
axes2 = axes.twinx()
axes2.plot(volume, salt, color='blue', label='Na+ (salt)', alpha=0.6)
axes2.set_ylabel('Concentration (mM)')
axes2.legend(loc='center left')
plt.show()