I was trying to use the generalized ion exchange (GIEX) model to simulate the elution curves under the pH effect for mixed-mode chromatography. Before I moved to parameter estimation, I used GIEX to create the reference model to obtain the h5 file. The experiment was done by step elution with pH values from 4.6 to 6.6. The first component of the GIEX model is salt. Since the salt concentration does not change during the experiment operations, the parameters for salt are set to 0. As I don’t have prior knowledge about the GIEX model parameters, I just set up some parameters to make sure the code is working, but the error attached below pops out.
I also uploaded the code to run the GIEX model. I am not sure whether the failure comes from the wrong parameters or the model setup. Would you please take a look at the code and let me know where the error comes from?
Glad to hear of another person working on modeling for multimodal chromatography! I recently published a paper showing the GIEX model is appropriate for modeling both MMCEX and MMAEX (see pre-proof Isotherm model discrimination for multimodal chromatography using mechanistic models derived from high-throughput batch isotherm data - ScienceDirect). If you want some examples on the model parameters, you can take a look at the word document in the SI section which has the parameters for both Capto MMC and Adhere using a variety of isotherm formalisms, where SMA Ext. is the same as GIEX. Keep in mind that we use a reference pH (pH = pH_actual - pH_reference) of 5.5. This helps with scaling of the pH-dependent parameters and also makes the changes w.r.t. pH more obvious. For example, at pH_actual = pH_reference the effective Keq and nu are equal to their base values (Keq0 and nu0). Depending on the value of the reference you use, it may result in negative pH values which may produce errors in the parameter value checks (not sure if this is present in CADET-Python, but we had to remove these checks in the CADET-MATLAB). All the simulations in this paper were done using ChromX, but we were able to reproduce the simulation results pretty closely in CADET—so they should be valid to use. Lastly, take care to mind the units of Ks and Kp! They are listed as M^-1 in the paper but need to be in m^3/mol (mM^-1) in CADET, so you will need to divide them by 1000. One of the findings in our paper was the importance of the batch data loading conditions—where, in short, the low concentration data should not be used when the initial slopes are high. Because of this, you should refer to the parameter values for Capto MMC BNL (batch no linear) and Capto Adhere BFI (batch full isotherm) using the SMA Ext. isotherm. Overall, this should give you a good idea of the expected ranges for isotherm parameter values for multimodal chrom.
I was also wondering if by pH effect you are referring to a general pH dependence in multimodal or the induced pH gradients (pH transients) seen for weak cationic groups. If you are referring to the later, I am also working on a paper w.r.t. to modeling the pH transients and the effects on protein elution using an in-house MATLAB script (adapted from Pabst and Carta 2007) coupled with CADET. I can share this code at some point (this summer once my PhD thesis is published). If you are referring to the general pH dependence, I think the GIEX model is pretty effective in capturing this—although for weakly charged ligands the transition to a very hydrophobically-dominated regime (due to the resin becoming less charged) and the presence of pH transients makes modeling the system quite difficult. From my studies so far, you cannot obtain accurate step elution predictions without accounting for the pH transients (assuming they are present in the system). They result in a forward shifted elution peak in the salt step which is impossible to predict assuming a constant pH. Worth noting that you can abrogate the pH shift by using an appropriate buffer at a high concentration (see Pabst and Carta 2007 in JChromA).
As for the script you shared… I am not familiar with this error coming from your simulation. From a quick glance, the model setup looks okay, but there are a number of errors in the specified inputs. Here is a list below. To obtain a general idea of the isotherm parameter values, you can refer my paper linked earlier in this post.
Particle diffusion for pH and salt are way too small. Should be around 100x higher than that of the protein.
Interstitial velocity is specified in terms of cm/s and should be in m/s
I generally don’t use the kD parameters since it is convenient to set kD = 1 such that kA is equal to Keq. For the kD lin, quad, salt, protein you can keep these all set to zero.
Lambda (ionic capacity) is much too small. It is in terms of mol ions / m^3 solid phase. One example calculation of this: Capto MMC having Ee = 0.4 and Ep = 0.91 (Et = Ee + (1-Ee)*Ep = 0.946) with ionic capacity value in the middle of the range reported from Cytiva (80 mol / m^3 total column volume). Here you would take 80 / (1-0.946) = 1481 mol/m^3 solid phase to be the value for CADET. This value is heavily dependent on the porosity values you have and is critical to get right, especially if you ever plan to use another resin lot in your studies.
Not sure how you decided on those reference concentrations, but I would keep them equal to 1 for now until the other problems are fixed. This is optional though.
Initial c and q for salt are set to zero. This will break the ion exchange model because having zero salt results in dividing by zero in the isotherm which might be the reason for your error, but I am not sure about that. The initial c for salt should be the concentration of salt in your equilibration buffer and the initial q is lambda (assuming no protein is bound).
I think section continuity should be [0, 0].
Constant coefficients for salt are zero, refer to point #6. Protein concentration seems to be way too high. Did you convert it to units of mol/m^3 (1:1 with mM)?
The flow rate seems to be impossible given your cross sectional area and velocity. I calculated 4.4047e-8 m^3/s.
I think that is it but I could have missed something. In general, there are a lot of inputs that need to be specified and it is critical to be very careful in assigning these. It is especially important to always check your units because this is one of the most common mistakes that can make the simulation invalid even if you make one mistake. Take your time and review the inputs carefully!!
Let me know if the simulation runs after you make all these changes.
Thanks so much for your detailed explanation! The experiment for mixed-mode chromatography is conducted under pH gradient elution and step elution. So by pH effect, I am referring to the latter, where the protein is eluted under varying pH values. So in my model setup, the pH value keeps constant at 4.6 during the loading and equilibrate phase, followed by a step elution with pH increasing from 4.6 to 6.6.
Thanks for your input about the model setup. As I am a new to IEX or MMC, I am not familiar with the parameter values. I will revise the model based on your suggestion and let you know if the simulation works or not.
Your help is much appreciated! Thanks again for your help!
Thanks again for your input. I have revised the errors based on your suggestions. In order to get the code running, I used the parameters given in your paper (SMA Ext. BFI and BNL) to simulate the process, but I got the same error as before. Would you give me some hints about how to get the code running?
Not sure if this will fix the error (I have never seen this error before), but I noticed that your constant coefficient terms still have salt concentration equal to zero. This should be changed to reflect the counterion concentration in the buffer, even if you are doing pH-based elution. Having zero salt causes problems in SDM-type models because it results in a divide by zero error (see Cs^-nu term).
After a quick glance, it seems like the rest of the code is correct.
the error arises from the fact that you specified the initial condition of the stationary phase for the second component of the Generalized Ion Exchange. The second component represents a non-binding modifier, hence no initial condition is required. If you specify the initial condition of the stationary phase as follows, the simulation should not fail at least for me the simulation runs.
