Mechanistic modeling, simulation, and optimization of mixed‐mode chromatography for an antibody polishing step

Mixed‐mode chromatography combines features of ion‐exchange chromatography and hydrophobic interaction chromatography and is increasingly used in antibody purification. As a replacement for flow‐through operations on traditional unmixed resins or as a pH‐controlled bind‐and‐elute step, the use of both interaction modes promises a better removal of product‐specific impurities. However, the combination of the functionalities makes industrial process development significantly more complex, in particular the identification of the often small elution window that delivers the desired selectivity. Mechanistic modeling has proven that even difficult separation problems can be solved in a computer‐optimized manner once the process dynamics have been modeled. The adsorption models described in the literature are also very complex, which makes model calibration difficult. In this work, we approach this problem with a newly constructed model that describes the adsorber saturation with the help of the surface coverage function of the colloidal particle adsorption model for ion‐exchange chromatography. In a case study, a model for a pH‐controlled antibody polishing step was created from six experiments. The behavior of fragments, aggregates, and host cell proteins was described with the help of offline analysis. After in silico optimization, a validation experiment confirmed an improved process performance in comparison to the historical process set point. In addition to these good results, the work also shows that the high dynamics of mixed‐mode chromatography can produce unexpected results if process parameters deviate too far from tried and tested conditions.


| INTRODUCTION
The purification processes for monoclonal antibody production in the biopharmaceutical industry commonly consist of a series of chromatography steps, starting with product capture using Protein A chromatography, followed by two or three polishing steps. 1 In the past, ion-exchange chromatography (IEC) or hydrophobic interaction chromatography (HIC) was mainly used as a polishing step, but in the last decade, mixed-mode chromatography (MMC) has been increasingly considered as an alternative. 2 Mixed-mode resins offer multiple binding interactions and have been demonstrated to be suitable for the purification of antibodies. 3,4 The MMC ligands often include features of IEC and HIC ligands. By adjusting the operational conditions, the binding modes can be reduced to a single interaction 5,6 or used together to affect the separation. MMC is highly suitable for application in pH-tunable hydrophobic 7 or high salt conditions, 8 and for capturing specific proteins under physiological conditions. 9 The fundamentals of MMC were the subject of experimental studies, 10 as well as simulative approaches 11 that generate valuable insights into protein-ligand interactions, but also showed that the inherent complexity necessitates further research. Despite the open questions, macroscopic models have been successfully used for process modeling. The Steric Mass Action model (SMA) 12 originally developed for IEC has been employed successfully for MMC. 13 Based on SMA and Mollerup's thermodynamic framework, 14,15 Nfor et al. 16 developed a model for protein adsorption in IEC/HIC MMC promising broader applicability. Further studies extended this model to describe observed adsorption behavior. Salt-dependent pore diffusion was employed for high-salt binding 17 and changes in water structure 18 have been taken into account, similar to the usage in pure HIC models. 19 A recent study performing model discrimination 20 has shown that a high level of model complexity is not necessary in every case and that, for example, the water structure change can be neglected. However, the same study also showed discrepancies of the model by Nfor  of the surface saturation function is independent of the chromatography mode and is thus assumed to be transferrable to MMC. At the same time, the model complexity is reduced, which reduces the risk of overfitting with a small data set for model calibration.
This study showcases application of our new modeling approach to the optimization of a MMC step in the purification process for an industrial antibody. Through our case study, we aim to contribute to the evaluation of the methodology's capabilities and limits, since this application has been considered a challenge in recent industrial position papers. 23,24 Following, References 25,26 the experimental design for model calibration is chosen depending on desired model capability. In this study, a small experimental data set is used for process optimization, which shall predict trends sufficiently well to enable the identification of an improved process set point. This differs from the procedure of using simulations as a substitute for real experiments in process characterization studies. In case a model with greater predictive power is to be created, the experimental design typically includes linear gradient elutions at low load density and possibly experiments with spikedin impurities 27,28 to put more emphasis on accurately describing the trace impurities in particular. Other studies used experiments under non-binding conditions or breakthrough curves to describe the fluid dynamics of the system in more detail. 29 All of these additional experiments improve the prediction quality of the model but were not con-  16 developed a model for protein adsorption in IEC/HIC MMC by assuming that both adsorption modes occur simultaneously. A mole of protein molecule P in solution is assumed to bind to ν moles of IEC ligands SL and n moles of HIC ligands L, forming 1 mole of proteinligand complexes PL, and exchanging ν moles of salt counter-ions S: Applying the law of mass action, the equilibrium formulation was derived to be where, c and c salt are the molar protein and salt concentration in solution, respectively, and q is the molar protein concentration bound to the stationary phase. The molarity of the solution in the pore volume c υ is assumed constant and the protein solute activity coefficient is modeled as with constant interaction parameters K S and K P .
In the linear isotherm range, the normalized concentration of available IEC ligands q IEC and HIC ligands q HIC can both be assumed to be one and K P can be set to be zero, as its influence is negligible.
Including c υ in the equilibrium constant k eq , we obtain the following isotherm function F eq for the linear range: The exponentiated capacity Λ ν was hidden in k eq after normalizing the available ligand concentration. It is factored out as ν is considered pH dependent. If the model was to be used for different resin lots with varying capacities, the exponentiated HIC capacity Λ n remaining in k eq would need to be factored out as well.
Following Reference 31 for IEC and Reference 32 for HIC, the pH dependency of the equilibrium constant is empirically included as Further, according to References 22,31, the protein charge ν's dependency on the pH is included as a second-order polynomial around a reference value pH ref .
The general structure of the non-linear isotherm follows the CPA formalism for IEC. 21 The rate of change of the bound protein concentration is described by a constant kinetic rate k kin and adsorption and desorption terms, multiplied by the protein concentrations in solution and adsorbed state, respectively.
The adsorption term consists of the equilibrium function defined above and the available surface function B from the CPA model, which depends on the specific surface area A s and colloid radius a and is explained in detail in Reference 22.

| Column model
In line with other industrial case studies, 26,30 the transport dispersive model (TDM) 33 was selected to describe the temporal change of the solute bulk concentration c i of solute i in the interstitial column volume where, x represents the axial position within the column, t is the time, D ax denotes the axial dispersion coefficient, k eff,i is the effective mass transfer coefficient, u int is the interstitial velocity, ε represents the void fraction, and c p,i is the concentration of the ith solute in the pore volume. As this model only considers an average concentration in the pore volume, k eff,i accounts for both internal and external mass transfer resistance. A General Rate Model (GRM) 33 can be used to separate these two mass transfer effects and to achieve better predictions, especially under overloaded conditions. In order to determine the additional parameters of the GRM, however, dedicated experiments are necessary, ideally breakthrough curves, 29 which were considered out-of-scope for this study.
The mass transfer equation for the pore volume is given by where, ε p denotes the particle porosity, d p is the adsorber bead diameter, and q i represents the concentration of the ith solute with respect to the adsorber skeleton volume.The column model is complemented with Danckwerts boundary conditions where, c in,i is the prescribed concentration of species i at the inlet of the column.

| Chromatographic instrumentation
Experiments were performed on an ÄKTA

| Buffers and feedstock
To operate the chromatography in pH-controlled bind-and-elute mode, a 25 mM sodium phosphate buffer at pH 7.7 was used for equilibration and wash. Elution was performed with a 50 mM sodium citrate buffer with an additional 200 mM NaCl at pH 3.4. The different pH profiles were mixed from these two buffers.
The monoclonal antibody used in this study is derived from industrial CHO cultivation and is of IgG1 subclass. The sample used is a Protein A eluate after virus inactivation titrated to pH 6. It includes a monomer species, high molecular weight (HMW) and low molecular weight (LMW) variants, and host cell proteins (HCPs). The HMW and HCP species were considered impurities to be removed in this case study.

| Experimental design
Historical data showed that an intermediate NaCl concentration, at the minimum of IEC and HIC interactions, alone is not sufficient for elution. However, a pH step to 4.2 had historically shown good results (cf. Table 7). Based on this information, six experiments were planned for model calibration. The first run aims to evaluate the behavior at a potential process set point in which the column is loaded to 85 g/L resin load density and then eluted over 5 column volumes (CV) in a pH 4.2 step. Runs 2, 3, and 4 instead elute in a 10, 20, and 30 CV pH gradient and are meant to determine the change of elution pH as a function of gradient slope, similar to experimental designs used for Yamamoto's method. 36

| Numerical methods
The simulations were performed using the software package ChromX (GoSilico GmbH, Karlsruhe, Germany). 38  Step 1 + step 2 70 - Step 1: 4.5 Step 2: 4.2 Step 1: 0.32 Step 2: 0.29 19 stepping algorithm. 39 The resulting linear systems are solved with UMFPACK. 40 pH was simulated as a mobile phase modifier according to Equation (8). To account for the non-linearity of the pH gradient, a transformation of the inlet concentration c in in Equation (10) was applied for the pH component. A fourth-order polynomial was used to match the pH simulation to the measured pH signal.

| Model parameter estimation
The model parameter estimation process was performed sequentially.
First, the dominating monomer species was modeled as a single-

| Model-based process optimization
After model calibration and quality assessment, the model was used for in silico process optimization. The process parameters considered for optimization were the load and elution conditions, the flow rate, and eluate pooling conditions. The chosen parameter ranges are shown in Table 5. The boundaries for load density and flow rate include process options outside the calibration space to explore further scenarios.
Adaptive simulated annealing was used to systematically vary the selected process parameters within defined boundaries and minimize the compound objective function shown in Equation (12).
The first two terms aim to maximize product yield and purity, with purity being weighted by a factor of 10. The third term penalizes a high HMW content. The fourth and fifth term reward higher recovery and thereby load density, as well as a small pool volume.
From the resulting process options, two conditions were selected to be confirmed experimentally.
The first selected run used the same feed material as the calibration runs. Optimization resulted in the same elution step at pH 4.2 but at a lower elution conductivity. While the process flow rate remained almost unchanged, the load density was reduced compared to that in calibration run 1. These conditions promised an increased product yield at improved HMW clearance.
Further, a second optimized process condition was selected where HMW depletion was favored over higher product yield. An

| Sample properties
The feed analysis by SEC showed 95.9% monomer, 3.3% species with lower molecular weight, and 0.8% species with higher molecular weight. After careful consideration of both the SEC data and the To calculate the molar concentrations of the species and approximate their colloid radii, the monomer's molecular weight was assumed to be 150 kDa. Accordingly, the LMW and HMW molecular weights were assumed to be 75 and 300 kDa, respectively. A colloid radius of 5.5 nm was used for the monomer; for the LMWs and HMWs, the radii 4.1 and 7.2 nm, respectively, were estimated using the empirical relation given in Reference 22. The extinction coefficient of the monomer was known, and the values for HMW and LMW were estimated from the peak areas of the pseudo-chromatograms generated from the offline analytics.

T A B L E 3 Summary of system and column parameters
It can be assumed that the HCPs are several species. However, because of their similar behavior, they were described as one component in the model. A molecular weight of 300 kDa and a radius of 7.2 nm were assumed for the HCPs, but this assumption has no noticeable effects on the model prediction due to the low HCP concentration in the feedstock (less than 0.001%).

| Model calibration
First, the pH transformation at the inlet was performed using a fourth-order polynomial to account for the non-linearity during gradient elution. As the step, elutions showed a rounded profile, a continuous stirred-tank reactor (CSTR) with a mixing rate of 0.0185 s À1 was added at the inlet to mimic the back-mixing effect of the buffer mixing chamber. Exemplary pH profiles are shown in Figure 1.
The remaining model parameters were determined by curve fitting using the measured UV signals at 280 and 300 nm and analytical fraction data of all runs following the stepwise approach outlined  Comparisons between fitted model curves and measured UV, as well as with fraction data are shown in Figures 2 and 3, respectively.
All model parameters are listed with their approximate confidence intervals in Table 4. The reference pH is in general arbitrary. In this case, it was set to 7.7, the pH of the washing buffer.
In general, the monomer species is described very well by the model. The monomer peak in the second elution step of run 6 is binding capacity is reached in these runs, which means that a more complex pore model would be necessary for better model performance. However, because later optimization was expected to determine a column load density which was below the dynamic binding capacity, the model complexity was not increased.

T A B L E 5 Process parameter ranges used for optimization and conditions chosen for validation experiments
The HMW species are well described with two components. The LMW and HCP species are not as well matched as the HMWs. An overprediction of both species is visible, but there are not enough data points to evaluate the cause of the deviation for those species.
For the purpose of process optimization, it is sufficient that the trends of the pool composition are described sufficiently well, which is the case for this model. For LMW and HCP, elution is rather triggered by the reduction in the equilibrium coefficient with decreasing pH. It is possible that a correlation between the equilibrium parameter and the charge persists that could not be resolved from the available data. Yamamoto runs as in Reference 35 would be necessary for improvement. The fact that HCP and LMW parameters are similar is because the several of the HCP data derived from the fractions were below the limit of quantitation for the assay and consequently has not resulted in an improvement in simulating the HCP separation behavior when optimizing parameters values algorithmically.
The given approximate 95% confidence intervals were determined in ChromX based on forward sensitivities and inversion of the Fisher information matrix. 17 The confidence estimates of the main component's parameters were determined using the L2 error norm.
As the impurities have little influence on the overall chromatogram, their sensitivities are low. For this reason, the confidence estimates of the impurity parameters were calculated based on the sensitivity with respect to the normalized residual mean square error. As shown in Table 4, the 95% confidence interval for certain parameters are still large. This does not indicate that the parameters are wrong or highly uncertain, but solely that a deviation of the parameter value has comparably little influence on the outcome of the simulation.

| Process optimization
The results of the two numerical process optimizations are shown in Table 5. The experiments were performed on the same column as that used for the calibration experiments. Two aspects were considered when comparing simulation and experimental results, the peak positions and shapes for all components, and the predicted versus experimental objective values. The visual comparison is shown in Figure 4.
Overall, the monomer predictions align well with the measurements.
The predicted tailing of validation run 1 is not as strong as observed in the measurement. This trend can also be seen for validation run 2, although the prediction is matching the measurement better for this experiment. The breakthrough behavior in validation run 2 is reflected correctly regarding the product amount that is breaking through, although the model is not able to describe the unusual shape of the breakthrough. The strip peaks are not covered by the monomer in the model.
The LMW and HMW species' elution peaks are predicted well apart from the HMWs in the strip peak. As also observed during model calibration, the model predicts a small breakthrough of HMW 1 which cannot be observed in the measurement (see calibration results in Figure 3).
The HCP elution profiles show that the HCP peaks are generally eluting with less tailing than predicted by the model. This was already observed during model calibration ( Figure 3). However, this deviation was not expected to strongly affect the model's ability to predict the HCP content of the elution pool, as the HCPs are eluting mainly with the product.
A comparison of the experimental eluate analysis and predicted pool content is given in Table 6. The data shows that the model was able to predict most objectives within expected accuracy. For validation run 1, the yield is predicted to be higher than experimentally confirmed. This is caused by the inaccurate prediction of the monomer tailing. As the model predicts less tailing and UV-based pooling is used, the recovered product in the pool is overpredicted. This also leads to an underprediction of the pool volume. The product purity is predicted correctly, as the contents of HMW, LMW and HCP are overall predicted well. HMW depletion is described particularly accurately.
For validation run 2, the step yield of all components is predicted especially well, which shows that the salt-dependent breakthrough behavior is captured by the model. Due to the slight difference between predicted and actual tailing behavior (Figure 4), the pool volume is underestimated by the model. The predicted positive effect of increased ionic strength during loading on HMW clearance could not be confirmed. Figure 4 shows that the model predicts the HMW 2 species to elute mainly in the strip, whereas HMW 1 is showing some breakthrough. This leads to the HMW amount in the pool to be underpredicted. The prediction of HCP clearance in validation run 2 and therefore content in the product pool is in good agreement despite visible differences in the peak shape in Figure 4h. The LMW content is slightly overpredicted by the model, as the LMW clearance in validation run 2 could be improved by the increased ionic strength in the feed. The yield of all components is also influenced by inaccuracies in the overall elution peak shape as the pooling is triggered based on the UV signal. This translates into deviations of predicted HMW and LMW yield of À26% and +22%, respectively.
It must be kept in mind that all calibration experiments were done without addition of NaCl to the feedstock and that the NaCl's influence was only visible during desorption. Therefore, validation run 2 was an extrapolation from only one data point that was not able to increase HMW clearance but only reduced the yield. In conclusion, it does not qualify as a prospective process set point. If validation run T A B L E 7 Comparison of historical experimental data and experimental results of validation run 1

Parameter
Historical data Validation run 1 Step

| CONCLUSIONS
The optimization of an antibody polishing step using a MMC resin  Table 7. It is to be noted that the historical runs used a column with a 20 cm bed height at the same linear flow rate as that utilized in these experiments. As a consequence of reducing the bed height to 10 cm but maintaining the same linear flow rate in these modeling experiments, the residence time was reduced by half, increasing the productivity of the process by a factor of two.
In summary, this case study showed that process modeling and