A General Rate Model of Ion-Exchange Chromatography for Investigating Ion-Exchange Behavior and Scale-up

This work presents a general rate model for ion-exchange in column chromatography with the steric massaction isotherm that is suitable for ion-exchange systems involving macromolecules such as large proteins. This comprehensive model considers axial dispersion, interfacial film mass transfer between the bulk-fluid phase and the particle phase, and intraparticle diffusion. The model system has been solved numerically using the finite element method and the orthogonal collocation method with a graphical user interface for Windows based personal computers. The effect of pH is modeled by including buffers, acids, and bases as species in the model system, allowing the induced pH gradients resulted from H+ or OHsorption by the ion exchanger to be described in the model. Feed profiles of all the species can be specified individually allowing complicated gradient patterns. The software, known as Chromulator-IEX, may be a useful tool for the investigation of ion-exchange chromatography behavior and its scale-up. A General Rate Model of Ion-Exchange Chromatography for Investigating Ion-Exchange Behavior and Scale-up


Introduction
In the biotechnology industry, downstream processing of a protein product usually takes up the majority of the overall production cost [1]. Because of a low product concentration in the feed and a high purity requirement, a multi-stage process is almost invariably required for a protein product intended for pharmaceutical applications. A typical downstream process consists of four stages, namely removal of insolubles, isolation, purification and polishing [2]. The core of the process is the purification stage. It centers on one or more liquid chromatography (LC) steps. The steps before these LC steps are used to disrupt the cells if the product is intracellular, to remove insolubles (e.g., cell debris), to eliminate the majority of the impurities and to reduce the feed volume in order to maximize the resolution and throughput of the LC steps while reducing costs by allowing the use of a smaller column volume. Among the various forms of LC, ion-exchange chromatography (IEC) is one of the most common because of its good resolution and high capacity. Reverse-phase liquid chromatography (RP-LC) typically achieves a better resolution than IEC. However, RP-LC uses a solvent mobile phase that can lead to denaturation of a protein product. Its loading capacity is also usually much lower than that of IEC.
To dislodge proteins from the ion-exchange resin inside an LC column after loading, a salt solution, often in a gradient feed pattern, can be used to exchange with the sorbed proteins. Another common practice is to use a pH gradient to elute the proteins with or without a salt gradient. Modeling IEC with pH gradients requires the concentration of H + (or OH -) to be known throughout the column. This is more complicated than modeling other ions because H + and OHcan combine with each other and with other ions in buffers in addition to binding to the resin. In order to calculate the pH, a charge balance must be performed on all ions in a particular location in the column. When buffers are present, a set of implicit equations must be solved for the H + or OHconcentration. If sorption of H + or OHis ignored, then the charge balance is not necessary, and the concentrations of coions (not including buffers) do not need to be known. However, induced pH transitions can occur as a result of this sorption. Ghose et al. [3] observed pH dips as large as 1 pH unit in the effluent of a cation-exchange column in response to step changes in NaCl feed, and a similar pH increase in anion exchange. In order to model these phenomena, sorption of H + and OHon the ion-exchange resin must be considered. A general rate model for LC consists of a bulk-fluid phase mass transfer partial differential equation (PDE) and a particle phase mass transfer equation coupled with an equation to describe the eluite-stationary binding mechanism (e.g., an isotherm). For high performance LC (HPLC), the concentration profiles for eluites in the bulk-fluid phase are very stiff. This makes solving the bulk-fluid phase PDE very demanding numerically. An efficient discretization method is needed for the column length axis. Raghavan and Ruthven used the orthogonal collocation (OC) method for fix-bed adsorption column modeling [4]. This method is efficient for relatively non-stiff systems. However, it cannot be used for stiff systems that require an excessive number of OC points because the OC method is designed to work stably for a few OC points. To overcome this limitation, Yu and Wang used the OC on finite element (OCFE) method [5] for ion-exchange chromatography modeling. Since there is no limitation on the number of finite element, the OCFE method can be used for very stiff concentration profiles. Wang's group subsequently used this method to simulate an IEC system with a mass-action isotherm in both equilibrium and nonequilibrium forms [6]. None of these models have included the calculation of pH in the presence of a buffer. Frey and coworkers [7] used a lumped pore resistance model to describe protein IEC with adsorbing buffers. This model was used to simulate the elution of bovine serum albumin (BSA) and hemoglobin on the Q-Sepharose FF (strong anion exchanger) gel in the presence of N-morpholinoethanesulfonic acid, acetate, and formate buffers [8].
Gu et al. investigated the size exclusion effect that could result in uneven (molar based) saturation capacities in a binary adsorption system with each species satisfying the Langmuir isotherm individually [9]. Their binary system had a large molecule that could not enter some small pores with binding sites inside and thus it had a smaller saturation capacity than the small molecule in the binary system. Such steric hindrance can also happen when a large molecule binds with an active site, but blocks the access to unoccupied active site(s) underneath it. They derived an extended Langmuir isotherm for uneven saturation capacities by introducing a discount factor. Their simulation work presented a "peak reversal" phenomenon as a result of relative selectivity reversal due to isotherm crossover in a binary adsorption chromatography system. For ion-exchange, Brooks and Cramer [10] devised a steric mass-action (SMA) isotherm to extend the mass-action isotherm that is commonly used for ion-exchange, by accounting for steric hindrance by sorbed macromolecules such as proteins that block the access to salt ions on the ion exchange resin by other macromolecules that want to replace the salt ions through ion-exchange. Their SMA derivations also lead to selectivity reversals, but not adsorption azeotropes in displacement chromatography. A theoretical model for the changes in SMA parameters of a protein in response to changes in the mobile phase pH was developed by Bosma and Wesselingh [11]. They found that for BSA binding on the Q-Sepharose FF, the characteristic charge could be considered constant and the log of the equilibrium constant was a nearly linear function of pH. Conversely, Strong and Frey [8] considered the characteristic charge of BSA to be a linear function of pH and the equilibrium constant to be constant under varying pH on the same ion exchanger.

Mass-Transfer Model
A general rate model was used to model the mass transfer in IEC. This model considers axial dispersion, interfacial film mass-transfer resistance, and intraparticle diffusion in the multicomponent LC system. The dimensionless equations for eluite species i are as listed below for the bulk-fluid phase and the particle phase, respectively [12].
The dimensionless boundary conditions are In Eq. (6), c fi (τ) is the dimensionless feed concentration of species i as a function of dimensionless time, and its value determines the type of chromatography (e.g., isocratic elution or gradient elution) being employed, or the sample injection shape (e.g., a pulse). Time zero is the time a species first enters the column inlet. The constants used to nondimensionalize the concentrations, C 0i , are arbitrary, but their values are usually chosen to be the largest value fed to the column. In the model system, Pe Li , η i and Bi i are dimensionless mass transfer parameters. The other dimensionless mass transfer parameter, ξ i , is calculated from 3Bi i •η i (1−ε b )/ε b in which ε b is the bed voidage. The estimation of these mass-transfer parameters from existing masstransfer correlations were recently described in detail by Gu et al. [12] who supplied a spreadsheet for their calculations as a supplement to their publication.

Adsorption Equilibrium Model
Reaction (7) shows the replacement of Na + initially sorbed on a cation-exchanger resin (sorbed phase or stationary phase) by Ca 2+ in the aqueous solution, For protein separations, a positively charged region of a protein will bind with the resin by replacing the sorbed Na + . For an anionexchanger, the initially sorbed ion is an anion such as Clwhich can be replaced by a protein using its negatively charged region. Eq. (8) is a generalized ion-exchange reaction for both cation exchange and anion exchange, depicting the exchange of species i in the solution with species j initially sorbed on the stationary phase. The reaction is assumed to be an elementary reaction. C denotes a concentration in the solution based on liquid volume while Q concentration in sorbed phase based on the particle skeleton volume (excluding macropores). Because the ion-exchange reaction happens in the particle phase rather than the bulk-fluid phase, C is given a subscript p to be consistent with the notations in Equations (1) and (2). The υ i and υ j values are the absolute values of the charges carried by species i and j, respectively. They also serve as stoichiometric coefficients to balance Eq. (8), Eq. (8) leads to the following mass action isotherm without using activity coefficients to describe the ion-exchange equilibrium, (10) in which N s is the number of species in the system. A Qi expression is needed from the isotherm to get q i in Eq. (2). The isotherm given by Eq. (9) and Eq. (10) is completely specified when N s −1 independent equilibrium constants, K ij , and the exchange capacity, Λ, are given. For convenience, species 1 is a salt ion, such as Na + or Cl -. It is used as the reference for all of the equilibrium constants. Using the relation below, all equilibrium constants can be calculated from the N s −1 independent equilibrium constants.
Unlike the Langmuir isotherm used in most LC models [12], In general there is no explicit analytical expression for Q i . Thus, the partial derivatives, ∂Q i /∂C pj , which are needed to solve the particle phase PDE after orthogonal collocation discretization [13], become implicit as well. This hampers the numerical solution to the IEC model. One way to avoid this difficulty is to step back to the isotherm precursor, i.e., the rate expression that models the ion-exchange reaction in Eq. (8) as shown below, where k ij is the forward reaction rate constant for the exchange of species i in solution with species j in the sorbed (stationary) phase as shown in Eq. (8), and k ji the backward reaction rate constant for the exchange of species j in the solution with species i in the sorbed phase. The rate constants can be related to the equilibrium constants by the following definition, The ion-exchange equilibrium isotherm expressed by Eq. (9) is a subset of Eq. (12). When the rate constants are sufficiently large, local equilibrium is achieved such that ∂Q i /∂t=0, Eq. (12) degenerates to Eq. (9). Eq. (10) is satisfied throughout the solution procedure if it is satisfied in the initial condition for the sorbed phase. A dimensionless form of Eq. (12) is where Da ij and Da ji are Damköhler numbers (dimensionless). The Damköhler numbers are related to the rate constants by For a system with N s species (including the ion that is used to equilibrate the column initially), N s (N s −1)/2 rate constants or Damköhler numbers can be chosen independently of the N s −1 equilibrium constants. In the nonequilibrium model used by Ernest et al. [6], k ij = 0 for all i, j ≠ 1. The kinetic expression given by Eq. (14) does not suffer from the drawback of species 1 having to be present for the system to approach equilibrium. It is not generally practical to determine the actual values of the rate constants, but by selecting large values for Da ij , equilibrium will be achieved, and the simulation results will not be sensitive to changes to these Da ij values as long as they are sufficiently large (e.g., 1000). Fixing arbitrary large values for Da ij leads to large values for k ij . Given K ij values, k ji values and subsequently Da ji values can be calculated from Eq. (13) and Eq. (15b), respectively. By doing so, the kinetic model in Eq. (14) reduces to equilibrium isotherm during simulation without the aforementioned implicit isotherm problem that hampers the numerical solution to the model. The downside is that the ODE system from the discretization of Eq.
(1) and Eq. (2) is not coupled with algebraic isotherm equations, but rather with additional ODEs from Eq. (14). This increases the size of the overall ODE system to be solved using an ODE solver.
In the SMA theory, species 1 is the salt ion species initially sorbed on the ion-exchange resin for the subsequent exchange with another ion (e.g., a protein). A macromolecule such as a protein molecule can exchange with the salt ion and covers some ion-exchange sites still sorbed with salt ions underneath. This steric hindrance effect makes these sorbed salt ions unavailable for ion-exchange with other proteins [10]. To derive the SMA isotherm, one replaces Q j in Eq. (9) and the right-hand side of Eq. (12) with by setting j=1. Q 1 is the concentration of salt ion sorbed on the ion-exchange resin available for ion exchange after considering the steric hindrance is calculated from the following equation in which σ is the steric factor for species i (excluding the salt, i.e., i≠1), Eq. (16) shows that 1 Q is smaller than Q 1 . This is because Q 1 is the concentration of salt ion sorbed on the ion-exchange resin available for ion exchange when there is no steric hindrance [10].

pH Model
By allowing both the equilibrium constant and the characteristic charge to vary with the pH, the elution of proteins with pH gradients can be modeled. In the numerical solution computer code, the charge is considered a polynomial in (pH−pH ref ). Because the results of Bosma and Wesselingh [11] suggest that the equilibrium constant is a more linear function of the pH in the logarithm, ln(pH) will be specified as a polynomial in (pH−pH ref ). pH ref , the reference pH, is arbitrary because the same function can be expressed for any pH ref by adjusting the coefficients of the polynomial.
To include pH effects in the model, the pH must be calculated at each point in the bulk-fluid and particle phases. The pH at any point and time in the column is affected by the concentrations of buffer, acid, and base, and the amount of H + and OHsorbed by the resin. The equations for the calculation of pH are shown below for cation exchange. For anion exchange, pOH is calculated by equations similar to these.

Numerical Solution
An analytical solution to a general rate model is only possible if the isotherm is a linear one as demonstrated by Lee et al. using Laplace transfer [14]. The model system consisting of Equations (1), (2) and (14) is highly nonlinear due to the complicated ion-exchange binding mechanism. In this work, it was solved using a numerical solution strategy described as follows. The finite element method is used to discretize the dimensionless column axis z and the OC method is used to discretize the dimensionless particle radius r [13]. For a typical IEC system, one or two interior OC points together with a built-in exterior OC point are sufficient to describe the concentration profiles in the particle phase. For large beads in some preparative-and large-scale applications, more OC points are needed as recently demonstrated by Gu et al. in adsorption LC [12]. After the discretization of Equations (1) and (2), the resultant ordinary differential equation (ODE) system will consist of N s (2N e +1)(2N c +1) ODEs if Ne number of finite elements and N c number of interior OC points are used for a system with N s species. The ODE system is solved using an efficient public-domain ODE solver known as VODE developed by Brown et al. [15]. The numerical solution in this work was implemented using the C++ computer language with a graphical user interface (GUI) to produce a software program known as Chromulator-IEX. The computation time i on today's personal computers is typically a matter of seconds. Very stiff systems may take up to several minutes or more (Figure 1). Figure 2 shows an example of the GUI of the software in which the absolute tolerance (atol) for the ODE solver is set to 1×10 -5 . The system has four species with species one being the salt ion that is used to elute out three other species isocratically. Species 1 is also used to equilibrate the column before injecting the feed containing species 2, 3 and 4. The pulse feed has duration of 2 in dimensionless time. It does not contain species 1. The results are shown in the simulated chromatogram in Figure 3 exhibiting the effluent concentration profiles of the four species. In Figure 2, the value of 1 in the c i (subscription i for initial rather than species i here) data entry column indicates that the IEC column is initially equilibrated with the C 0 concentration for species 1 (the salt ion), while the three zeroes point out that the IEC column does not contain any amount of species 2, 3 and 4 initially. If "Calculate pH" is checked in Figure 2, pKa and SF (salt fraction) values must be entered.

Results and Discussion
The software allows individual species' feed profiles to be specified separately. For example, the isocratic elution case in Figure 2 can be switched to a gradient elution with a salt (species 1) gradient from 0.5 dimensionless concentrations (i.e., half of the maximum in the feed) to 1 with a dimensionless time frame of 0 to 18. Figure 4A shows that "Gradient Elution, last component is displacer" and "Custom Mode of Operation" options are selected, which allow feed profiles to be specified in a new pop-up window. Figure 4B shows that species 1 concentration in the feed is 0.5 to 1 dimensionless concentration for the dimensionless time frame of 0 to 18. Figure 4C specifies that the feed for species 2 is a pulse injection lasting 2 dimensionless time. Species 3 and 4 are specified the same way as that for species 2. Figure  5 shows the simulated chromatogram using the parameters in Figure 4. The column is initially equilibrated with species 1 with a dimensionless concentration of one, which shows up in the effluent profile for the salt at dimensionless time 0. Compared with the isocratic elution outcome in Figure 3, the gradient elution in Figure 5 does not exhibit a drastic improvement, because the gradient turns out to have delayed the speed of peak migration causing the peaks to diffuse more. More elaborate feed patterns can be simulated by using the "Add" button Figure 4C to supply dimensionless concentration and dimensionless time data points piecemeal in order to construct a feed profile of a species.
To illustrate the transient ion-exchange behavior, the software utilizes a built-in movie player to show the bulk-fluid phase concentration profiles of all the species at different axial positions inside the IEC column. Figure 6 is a screenshot of the concentration profiles for the case in Figure 5 at dimensionless time τ = 4.5. The concentrations of the four species at column position z=1 (i.e., column exit) match those in Figure 5 at τ=4.5.
In Figure 7A, Ghose et al. [3] experimentally demonstrated the pH spike and pH dip phenomena by feeding three different buffers at the same pH to a Fractogel SO 3 -(M) ion-exchange column. The column was preconditioned with 2 M NaCl and 250 mM sodium citrate buffer at pH 5.5, equilibrated with 25 mM sodium citrate buffer at pH 5.5, and eluted at with 0.5 M NaCl and 25 mM sodium citrate buffer at pH 5.5. They commented that a pH dip below 5.0 or lower could cause significant bioactivity losses for some proteins. Figure 7B shows simulated results using Pe L =500, η=5, Bi=5 for all species, and K H+ , Na+ =1, Λ=95 equiv/L, ε b =0.4, ε p =0.5. The model in this work successfully demonstrates both the pH spike and the pH dip at the correct elution volume positions compared with the experimental data in Figure 7A. In a separate case, Figure 8 shows that the model fits the experimental data on Dowex 50 ion-exchange resin well. The literature data in Figure  8 are from the work by Dranoff and Lapidus [16].    A general rate model for IEC was presented and solved numerically. The model considered axial dispersion, interfacial film mass-transfer resistance, intraparticle diffusion and the steric mass-action isotherm for ion-exchange. The numerical difficulty caused by the implicit isotherm was avoided by using its precursor kinetic expression and setting the kinetic rate constants to sufficiently large values to achieve equilibrium. The model system coupled with the pH expression to calculate the pH value everywhere inside the column successfully simulated the phenomenon of induced pH transitions in IEC. Such transitions are important in chromatofocusing and in the prevention of protein denaturation due to large pH spikes and dips. The model and software presented in this work should be useful tools to those involved in the scale-up of ion-exchange chromatography. It can also be used in teaching and investigating IEC behavior.

Symbols
Bi = Biot number for mass transfer, kR p /(ε p D p )