Protein Concentration Not show up

CADET Troubleshooting
1

Hi,
I tried to simulate a two component (salt and protein) system. The initial protein concentration is 0.0115 mM in 20 mM salt. However, when I simulate inlet concentration, outlet concentration of CSTR (unit 2), inlet and outlet concentration of column, they all showed 0 protein concentration which isn’t reasonable. I’m wondering if I made any mistakes when setting up section and connections to cause this error. Here is my script to set up some basic values:

wash_start =    32                # s
grad_start =    wash_start + 365       # s
HighSalt_start = grad_start + 954       # s
t_end =         HighSalt_start + 500       # s

# Volumetric flow rate
vol_flow_rate = 1.6e-8                    # m3/s

# Salt inlet concentration 
salt_initial_concentration = 20         # mM
salt_final_concentration = 1020          # mM

# Lysozyme feed concentration 
protein_initial_concentration = 0.0115    # mM  

# Column dimensions (in meters)
column_dia = 0.007                          # m
bed_height = 0.025                        # m

lwe_model = get_cadet_template(n_units=5)

# Sections and Switches
lwe_model.root.input.solver.sections.nsec = 4
lwe_model.root.input.solver.sections.section_times = [0.0, wash_start, grad_start, HighSalt_start, t_end]
lwe_model.root.input.model.unit_000.sec_000.const_coeff = [salt_initial_concentration, 
                                                           0]
lwe_model.root.input.model.unit_000.sec_001.const_coeff = [salt_initial_concentration, 0]
lwe_model.root.input.model.unit_000.sec_002.const_coeff = [salt_initial_concentration, 0.0]
lwe_model.root.input.model.unit_000.sec_002.lin_coeff = [1.05,0.0]
lwe_model.root.input.model.unit_000.sec_003.const_coeff = [salt_final_concentration,0]

lwe_model.root.input.model.unit_001.sec_000.const_coeff = [salt_initial_concentration, 
                                                           protein_initial_concentration]
lwe_model.root.input.model.unit_001.sec_001.const_coeff = [0, 0.0]
lwe_model.root.input.model.unit_001.sec_002.const_coeff = [0, 0.0]
lwe_model.root.input.model.unit_001.sec_003.const_coeff = [0,0]

# Set the times that the simulator writes out data for
lwe_model.root.input.solver.user_solution_times = np.linspace(0, t_end, int(t_end) + 1)

# System of Unit Operations
n_comp = 2

#Feed
lwe_model.root.input.model.unit_000.unit_type = 'INLET'
lwe_model.root.input.model.unit_000.ncomp = 1
lwe_model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
#Inlet
lwe_model.root.input.model.unit_001.unit_type = 'INLET'
lwe_model.root.input.model.unit_001.ncomp = n_comp
lwe_model.root.input.model.unit_001.inlet_type = 'PIECEWISE_CUBIC_POLY'

#Mixer
lwe_model.root.input.model.unit_002.unit_type = 'CSTR'
lwe_model.root.input.model.unit_002.ncomp = n_comp
lwe_model.root.input.model.unit_002.init_volume = 5e-7
lwe_model.root.input.model.unit_002.init_c = n_comp*[0.0]

#Column Chromatogaphy Model
lwe_model.root.input.model.unit_003.unit_type = 'GENERAL_RATE_MODEL'
lwe_model.root.input.model.unit_003.ncomp = n_comp

# Column parameters
lwe_model.root.input.model.unit_003.col_length = bed_height
lwe_model.root.input.model.unit_003.cross_section_area = 3.142*((column_dia/2)**2)
lwe_model.root.input.model.unit_003.col_porosity = 0.268
lwe_model.root.input.model.unit_003.total_porosity = 0.873
lwe_model.root.input.model.unit_003.par_porosity = 0.827
lwe_model.root.input.model.unit_003.par_radius = 4.5e-5
lwe_model.root.input.model.unit_003.col_dispersion = 3.73e-7

lwe_model.root.input.model.unit_003.film_diffusion = [1.0E-5, 1.0E-5]
lwe_model.root.input.model.unit_003.par_diffusion = [6e-11, 8e-12]
lwe_model.root.input.model.unit_003.par_surfdiffusion = [0.0, 1e-10]

## OUTLET
lwe_model.root.input.model.unit_004.unit_type = 'OUTLET'
lwe_model.root.input.model.unit_004.ncomp = n_comp

## Discretization
set_discretization(lwe_model, n_bound=[0,1])

# Connections
lwe_model.root.input.model.connections.nswitches = 2
lwe_model.root.input.model.connections.switch_000.section = 0
lwe_model.root.input.model.connections.switch_000.connections = [0, 2, 0, 0, vol_flow_rate,
                                                                2, 3, -1, -1, vol_flow_rate,
                                                                3, 4, -1, -1, vol_flow_rate]

lwe_model.root.input.model.connections.switch_001.section = 1
lwe_model.root.input.model.connections.switch_001.connections = [1, 2, -1, -1, vol_flow_rate,
                                                                2, 3, -1, -1, vol_flow_rate,
                                                                3, 4, -1, -1, vol_flow_rate]

adsorption_model = "MULTI_COMPONENT_COLLOIDAL"
adsorption_parameters = Dict()

# Vijesh's model parameters for Bpp
adsorption_parameters.col_bpp_ph_exp = [0.0, 0.0] 
adsorption_parameters.col_bpp_salt_expargmult = [0.0, 0.0] 
adsorption_parameters.col_bpp_salt_expfact = [0.0, 0.2625] 
adsorption_parameters.col_bpp_salt_powexp = [0.0, -2.252] 
adsorption_parameters.col_bpp_salt_powfact = [0.0, 37220.0] 


adsorption_parameters.col_cordnum = 6.0 # hexagonal packing 
adsorption_parameters.col_kappa_const = 2.5 # nm
adsorption_parameters.col_kappa_exp = 0.0 
adsorption_parameters.col_kappa_fact = 0.0 
adsorption_parameters.col_kkin = [1.0e10, 1.0e10] # e6 to e12 for instant equilibration; 0.1-0.9 for kinetic effects -- doesn't work though
adsorption_parameters.col_linear_threshold = 1.0e-7 

# Vijesh's model parameters for Keq
adsorption_parameters.col_logkeq_ph_exp = [0.0, 0.0] 
adsorption_parameters.col_logkeq_salt_expargmult = [0.0, 0.0] 
adsorption_parameters.col_logkeq_salt_expfact = [0.0, -3.929] 
adsorption_parameters.col_logkeq_salt_powexp = [0.0, 1.715] 
adsorption_parameters.col_logkeq_salt_powfact = [0.0, 53430.0] 

adsorption_parameters.col_phi = 4.923e7 # m2 / m3
adsorption_parameters.col_radius = [1.59e-9, 4.5e-9] # make nonzero
adsorption_parameters.col_use_ph = False 
adsorption_parameters.is_kinetic = True 

lwe_model.root.input.model.unit_003.adsorption_model = adsorption_model
lwe_model.root.input.model.unit_003.adsorption = adsorption_parameters

lwe_model.root.input.model.unit_003.init_c = [0, 0]
lwe_model.root.input.model.unit_003.init_q = [0, 0.0]

run_simulation(lwe_model)
2

Hi,

I think I found the isssue.

First, you define dynamic changes in the connectivity.

# Connections
lwe_model.root.input.model.connections.nswitches = 2
lwe_model.root.input.model.connections.switch_000.section = 0
lwe_model.root.input.model.connections.switch_000.connections = [0, 2, 0, 0, vol_flow_rate,
                                                                2, 3, -1, -1, vol_flow_rate,
                                                                3, 4, -1, -1, vol_flow_rate]

lwe_model.root.input.model.connections.switch_001.section = 1
lwe_model.root.input.model.connections.switch_001.connections = [1, 2, -1, -1, vol_flow_rate,
                                                                2, 3, -1, -1, vol_flow_rate,
                                                                3, 4, -1, -1, vol_flow_rate]

But you also define dynamic changes to the concentration profile

# Sections and Switches
lwe_model.root.input.solver.sections.nsec = 4
lwe_model.root.input.solver.sections.section_times = [0.0, wash_start, grad_start, HighSalt_start, t_end]
lwe_model.root.input.model.unit_000.sec_000.const_coeff = [salt_initial_concentration, 0]
lwe_model.root.input.model.unit_000.sec_001.const_coeff = [salt_initial_concentration, 0]
lwe_model.root.input.model.unit_000.sec_002.const_coeff = [salt_initial_concentration, 0.0]
lwe_model.root.input.model.unit_000.sec_002.lin_coeff = [1.05,0.0]
lwe_model.root.input.model.unit_000.sec_003.const_coeff = [salt_final_concentration,0]

lwe_model.root.input.model.unit_001.sec_000.const_coeff = [salt_initial_concentration, protein_initial_concentration]
lwe_model.root.input.model.unit_001.sec_001.const_coeff = [0, 0.0]   #  <------- Look here
lwe_model.root.input.model.unit_001.sec_002.const_coeff = [0, 0.0]
lwe_model.root.input.model.unit_001.sec_003.const_coeff = [0,0]

The connectivity correctly switches to unit_001 this the concentration of that unit at that time is 0.

Once that’s fixed, and when you plot the solution, it still appears to be zero. But now you need to make sure that you adjust the axes accordingly since the salt concentration is so much larger than the protein. Took me a moment to find that “isssue” (which was no issue)… :sweat_smile:

Also, your example did not run correctly for me, there was an issue with consistent initialization. I just deactivated the isotherm, then it run.

Btw, I also prefer the approach of using multiple inlet units to describe the inlet profile of the column. Just wanted to let you know that if the concentration of your inlets stays constant, there is no need to explicitly set them for each section. It will just use the previous one.

So the following also works:

lwe_model.root.input.solver.sections.nsec = 4
lwe_model.root.input.solver.sections.section_times = [0.0, wash_start, grad_start, HighSalt_start, t_end]

lwe_model.root.input.model.unit_000.sec_000.const_coeff = [salt_initial_concentration, 0]

lwe_model.root.input.model.unit_001.sec_000.const_coeff = [salt_initial_concentration, protein_initial_concentration]
3 Likes
3

Hi Johannes,
Thank you so much for your solution and suggestion! It finally worked out!

Cheers,
Annie

4

Hi Johannes,
Sorry to bother you again. However, I realized a problem when I took a close look. If we did define section 001, 002, 003 for unit 001, when I graphed the inlet and outlet concentration profile of protein, it showed protein is continuously loaded, however, we only required to injected for 32s. In addition, I found out the salt concentration became 0 in the phase HightSaltPush, which we defined that the starting concentration will be “salt-final-concentration”, which should be 1020 mM.
Can I also please ask which script you used to deactivate isotherm?

Here is the updating script:

wash_start =    32                # s
grad_start =    wash_start + 365       # s
HighSalt_start = grad_start + 954       # s
t_end =         HighSalt_start + 1247        # s

# Volumetric flow rate
vol_flow_rate = 1.6e-8                    # m3/s

# Salt inlet concentration 
salt_initial_concentration = 20         # mM
salt_final_concentration = 1020          # mM

# Lysozyme feed concentration 
protein_initial_concentration =  0.0115   # mM  

# Column dimensions (in meters)
column_dia = 0.007                          # m
bed_height = 0.025                        # m

lwe_model = get_cadet_template(n_units=5)

# Sections and Switches
lwe_model.root.input.solver.sections.nsec = 4
lwe_model.root.input.solver.sections.section_times = [0.0, wash_start, grad_start, HighSalt_start, t_end]
lwe_model.root.input.model.unit_000.sec_000.const_coeff = [salt_initial_concentration, 0]
lwe_model.root.input.model.unit_000.sec_001.const_coeff = [salt_initial_concentration, 0]
lwe_model.root.input.model.unit_000.sec_002.const_coeff = [salt_initial_concentration, 0.0]
lwe_model.root.input.model.unit_000.sec_002.lin_coeff = [1.05,0.0]
lwe_model.root.input.model.unit_000.sec_003.const_coeff = [salt_final_concentration,0]

lwe_model.root.input.model.unit_001.sec_000.const_coeff = [salt_initial_concentration, protein_initial_concentration]

# Set the times that the simulator writes out data for
lwe_model.root.input.solver.user_solution_times = np.linspace(0, t_end, int(t_end) + 1)

# System of Unit Operations
n_comp = 2

#Feed
lwe_model.root.input.model.unit_000.unit_type = 'INLET'
lwe_model.root.input.model.unit_000.ncomp = 1
lwe_model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
#Inlet
lwe_model.root.input.model.unit_001.unit_type = 'INLET'
lwe_model.root.input.model.unit_001.ncomp = n_comp
lwe_model.root.input.model.unit_001.inlet_type = 'PIECEWISE_CUBIC_POLY'

#Mixer
lwe_model.root.input.model.unit_002.unit_type = 'CSTR'
lwe_model.root.input.model.unit_002.ncomp = n_comp
lwe_model.root.input.model.unit_002.init_volume = 5e-7
lwe_model.root.input.model.unit_002.init_c = n_comp*[0]

#Column Chromatogaphy Model
lwe_model.root.input.model.unit_003.unit_type = 'GENERAL_RATE_MODEL'
lwe_model.root.input.model.unit_003.ncomp = n_comp

# Column parameters
lwe_model.root.input.model.unit_003.col_length = bed_height
lwe_model.root.input.model.unit_003.cross_section_area = 3.142*((column_dia/2)**2)
lwe_model.root.input.model.unit_003.col_porosity = 0.268
lwe_model.root.input.model.unit_003.par_porosity = 0.827
lwe_model.root.input.model.unit_003.par_radius = 4.5e-5
lwe_model.root.input.model.unit_003.col_dispersion = 3.73e-7

lwe_model.root.input.model.unit_003.film_diffusion = [1.0E-5, 1.0E-5]
lwe_model.root.input.model.unit_003.par_diffusion = [6e-11, 8e-12]
lwe_model.root.input.model.unit_003.par_surfdiffusion = [0.0, 1e-10]
## OUTLET
lwe_model.root.input.model.unit_004.unit_type = 'OUTLET'
lwe_model.root.input.model.unit_004.ncomp = n_comp

## Discretization
set_discretization(lwe_model, n_bound=[0,1])

# Connections
lwe_model.root.input.model.connections.nswitches = 2
lwe_model.root.input.model.connections.switch_000.section = 0
lwe_model.root.input.model.connections.switch_000.connections = [0, 2, 0, 0, vol_flow_rate,
                                                                2, 3, -1, -1, vol_flow_rate,
                                                                3, 4, -1, -1, vol_flow_rate]

lwe_model.root.input.model.connections.switch_001.section = 1
lwe_model.root.input.model.connections.switch_001.connections = [1, 2, -1, -1, vol_flow_rate,
                                                                2, 3, -1, -1, vol_flow_rate,
                                                                3, 4, -1, -1, vol_flow_rate]

adsorption_model = "MULTI_COMPONENT_COLLOIDAL"
adsorption_parameters = Dict()

# Vijesh's model parameters for Bpp
adsorption_parameters.col_bpp_ph_exp = [0.0, 0.0] 
adsorption_parameters.col_bpp_salt_expargmult = [0.0, 0.0] 
adsorption_parameters.col_bpp_salt_expfact = [0.0, 0.2625] 
adsorption_parameters.col_bpp_salt_powexp = [0.0, -2.252] 
adsorption_parameters.col_bpp_salt_powfact = [0.0, 37220.0] 


adsorption_parameters.col_cordnum = 6.0 # hexagonal packing 
adsorption_parameters.col_kappa_const = 2.5 # nm
adsorption_parameters.col_kappa_exp = 0.0 
adsorption_parameters.col_kappa_fact = 0.0 
adsorption_parameters.col_kkin = [1.0e10, 1.0e10] # e6 to e12 for instant equilibration; 0.1-0.9 for kinetic effects -- doesn't work though
adsorption_parameters.col_linear_threshold = 1.0e-7 

# Vijesh's model parameters for Keq
adsorption_parameters.col_logkeq_ph_exp = [0.0, 0.0] 
adsorption_parameters.col_logkeq_salt_expargmult = [0.0, 0.0] 
adsorption_parameters.col_logkeq_salt_expfact = [0.0, -3.929] 
adsorption_parameters.col_logkeq_salt_powexp = [0.0, 1.715] 
adsorption_parameters.col_logkeq_salt_powfact = [0.0, 53430.0] 

adsorption_parameters.col_phi = 4.923e7 # m2 / m3
adsorption_parameters.col_radius = [1.59e-9, 4.5e-9] # make nonzero
adsorption_parameters.col_use_ph = False 
adsorption_parameters.is_kinetic = True 

lwe_model.root.input.model.unit_003.adsorption_model = adsorption_model
lwe_model.root.input.model.unit_003.adsorption = adsorption_parameters

lwe_model.root.input.model.unit_003.init_c = [salt_initial_concentration, 0]
lwe_model.root.input.model.unit_003.init_q = [0, 0.0]

run_simulation(lwe_model)

# Plot the simulation
fig, ax1 = plt.subplots()

# Plot the salt trace
ax1.plot(lwe_model.root.output.solution.solution_times,
    lwe_model.root.output.solution.unit_003.solution_outlet[:,0], 'k', label='Salt')
ax1.set_xlabel('Time(s)')
ax1.set_ylabel('Salt Concentration (mM)')

# Create second axis
ax2 = ax1.twinx()

# Plot the protein (mAb) trace
ax2.plot(lwe_model.root.output.solution.solution_times,
    lwe_model.root.output.solution.unit_003.solution_outlet[:,1], 'r', label='Lysozyme')
ax2.set_ylabel('Protein Concentration (mM)')
"""If you want to set your own y-axis limits for the protein trace in the
figure, un-comment the line below and set lower-limit, upper-limit respectively
"""
#ax2.set_ylim(0.00,0.6)

# Formatting
lines, labels = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2)

# Save a .csv file in the working directory (where this .py file is kept)
col_1 = lwe_model.root.output.solution.solution_times
col_2 = lwe_model.root.output.solution.unit_003.solution_outlet[:,0]
col_3 = lwe_model.root.output.solution.unit_003.solution_outlet[:,1]

output_df = pd.DataFrame({'time': col_1, 'salt_conc':col_2, 'protein_conc':col_3})
output_df.to_csv('output3.csv', index=False)

Thank you very much for your help!

5

Hi Annie,

ok, maybe my description wasn’t complete. :sweat_smile:

When the number of switches is lower than the number of sections, they will cycle.
E.g.

  • Section 0: Switch 0
  • Section 1: Switch 1
  • Section 2: Switch 0
  • Section 3: Switch 1

etc.

So I assume that might cause an issue. Following the Zen of Python, “Explicit is better than implicit” I would suggest always defining the corresponding switches for every section unless it’s trivial (i.e. they are constant for all sections).

So maybe this fixes your problem.

Also, to deactivate the isotherm, I just omitted the adsorption_model and the adsorption branch of the config.

Also Also, Maybe you need to adjust your initial conditions. I don’t know the colloidal isotherm but often when salts are involved, make sure that init_q is not zero (e.g. set the salt loading to the capacity of the column) and maybe even init_c (if you equilibrate with salt buffer). Otherwise consistent initialization has issues (which seems to be the case for your simulation).

To make our life a bit easier when providing help, I would like to ask you to format your code using code fences (three backticks (`) on a single line).

For example,

```
print(‘hello annie’)
```

will be rendered as

print('hello annie')

Moreover (and more importantly), please make sure that the scripts you provide are self contained. If variables are missing or functions are imported that we don’t have access too, it’s difficult for us.
To check, just run them in a fresh kernel and see if any errors come up (except the ones you want to tell us about, of course! :nerd_face:)

Let me know if you have more questions.

6

Hi Johannes,
Thank you very much for your detailed explanation! I will keep in mind in later posted questions.

Kind regards,
Annie