Prediction tool for loading, isocratic elution, gradient elution and scaling up of ion exchange chromatography of proteins

Accepted Manuscript Title: Prediction tool for loading, isocratic elution, gradient elution and scaling up of ion exchange chromatography of proteins Authors: Wojciech Kazimierz Marek, Dominik Sauer, Astrid Durauer, ¨ Alois Jungbauer, Wojciech Pi˛atkowski, Dorota Antos PII: DOI: Reference:

S0021-9673(18)30812-4 https://doi.org/10.1016/j.chroma.2018.06.057 CHROMA 359501

To appear in:

Journal of Chromatography A

Received date: Revised date: Accepted date:

6-3-2018 20-6-2018 22-6-2018

Please cite this article as: Marek WK, Sauer D, Durauer ¨ A, Jungbauer A, Pi˛atkowski W, Antos D, Prediction tool for loading, isocratic elution, gradient elution and scaling up of ion exchange chromatography of proteins, Journal of Chromatography A (2018), https://doi.org/10.1016/j.chroma.2018.06.057 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Prediction tool for loading, isocratic elution, gradient elution and scaling up of ion exchange chromatography of proteins Wojciech Kazimierz Mareka,*, Dominik Sauerb, Astrid Dürauerb,c, Alois Jungbauerb,c, Wojciech Piątkowskia, Dorota Antosa of Chemical and Process Engineering, Powstańców Warszawy Ave. 6, 35-959 Rzeszów, Poland bAustrian Centre of Industrial Biotechnology, Muthgasse 11, 1190 Vienna, Austria cDepartment of Biotechnology, University of Natural Resources and Life Sciences Vienna, Muthgasse 18, 1190 Vienna, Austria

*Corresponding

SC R

IP T

aDepartment

author; E-mail address: [email protected] (W.K. Marek), tel.: +48 177432020; fax: +48

U

178543655

Highlights

Mathematical tool for design and scaling up of protein chromatography is

N



A

suggested

Moment analysis was used for determining kinetic parameters of mass transfer



Two model ion-exchange chromatography systems were selected for the study



The design procedure was successfully verified over a wide space of the operating

CC E

Abstract

PT

conditions

ED

M



An efficient mathematical tool for the design and scaling up of protein chromatography is suggested, in which the model parameters can be determined quickly over a wide operating space without large material investments. The design method is based on

A

mathematical modelling of column dynamics and moment analysis. The accuracy of the dynamic models that are most frequently used for simulations of chromatographic processes is analyzed, and possible errors that can be generated using the moment analysis are indicated. The so-called transport dispersive model was eventually employed for the process simulations. The model was modified to account for the protein dispersion in void volumes of chromatographic systems. The manner of the model calibration was 1

suggested, which was based on a few chromatographic runs and verified over a wide space of the operating parameters, including composition and flow rate of the mobile phase, column dimensions, residence time, and mass loading. The model system for the study was ion-exchange chromatography. The analysis was performed based on the elution profiles of basic fibroblast growth factor 2 and lysozyme, on two different IEX media.

IP T

Keywords: mathematical modeling; protein chromatography; scaling up; fibroblast

SC R

growth factor; lysozyme

1. Introduction

U

The realization of chromatographic processes on a large production scale requires a design stage which involves decisions on the choice of operating parameters, such as:

N

column dimensions, composition and flow rate of the mobile phase, and feed load, with

A

the goal of defining satisfactory process performance, i.e., purity and cost of the operation

M

[1-3]. Ion exchange chromatography of proteins is very sensitive to the composition of salt in the feed stock and in the mobile phase when loading, washing and eluting columns.

ED

This influences the performance parameters such as dynamic binding capacity, elution volume, protein concentration in the elution peak, and peak retention [1,2,4]. To accurately describe and predict ion-exchange chromatography of proteins, the relevant

PT

parameters must be estimated as a function of salt and protein concentration. That goal must be however achieved rapidly and at low costs of material resources. Recently, high

CC E

throughput experiments methods (HTPE) are used, in which the process operating window is defined based on statistical analysis that exploits data acquired from microscale experiments using chromatographic columns with very small volumes, i.e.,

A

maximum 1 mL [5-7]. Nevertheless, the scalability of such experiments has some limitations. More satisfactory results can be achieved using the model-based design approach (MBD) that exploits mechanistic models of the process [8-11] or by efficient coupling of both HTPE and MBD methods [12]. Dynamic models that are commonly used in MBD to predict chromatographic processes account for axial dispersion, external and internal mass transport resistances, and kinetic rates of adsorption-desorption [13-17]. In heterogeneous models, all partial 2

processes are described directly by adequate differential mass balance equations for the liquid phase and the adsorbent particles, whereas in pseudo-homogeneous models, apparent parameters are employed which lump together contributions of different kinetic effects in mass balance equations for the liquid phase [18-20]. The models formulated in this way may be applied to several tasks: prediction and characterization of the chromatographic step in terms of the process performance, the selection of optimal operating conditions [16]; definition of critical process parameters and addressing of issues of scale-up and column scheduling [22].

IP T

[9]; the sensitivity and robustness analysis for quantifying the risk of batch failure [21];

SC R

The MBD method requires the determination of a number of kinetic and thermodynamic parameters underlying the dynamic model. Increased model complexity

results in an increase in the number of parameters to be quantified. Since the contributions of mass transport resistance and dispersion in the chromatographic column

U

to band broadening are additive, all kinetic parameters of heterogeneous models cannot

N

be determined independently using chromatographic elution experiments. Simple elution

A

studies often result in misinterpretation of the physical meaning of the model parameters and their dependencies on the operating conditions, and therefore, reducing the

M

operating space in which appropriate accuracy of predictions can be achieved. Moreover, the solution of complex heterogeneous models involves high computational efforts, which

ED

hinder their use for the process optimization or online control, where a fast-simulation software is required. On the other hand, model oversimplifications or incorrect use of

PT

apparent kinetic coefficients can result in erroneous predictions of chromatographic band profiles and failure in scaling up the process.

CC E

The accuracy of predictions of a chromatographic process by pseudo-homogenous models has already been examined based on the moment analysis method [18,23-25]. However, the analysis presented by these authors concerned mainly HPLC conditions; it did not cover the range of large values of the plate height, i.e., 0.1 cm and higher, which is

A

specific for protein separations in preparative and industrial scale applications [1]. Another issue, which is often neglected while modelling and scaling up

chromatographic processes, is the contribution of dispersion in void volumes of chromatographic system to the column efficiency. In case of mini-columns that contribution is significant, and it can determine the pattern of protein retention.

3

Therefore, it has to be distinguished from other kinetic effects accompanying protein elution. The goal of the present work was to develop a simple but efficient mathematical tool for the design and scaling up of protein chromatography, in which the contribution of void volumes to the process efficiency was accounted for. In the first part of the study, we present an analysis of the accuracy of the dynamic models that are most commonly used for simulations of chromatography processes of low and very low efficiency, including a

IP T

heterogeneous model, i.e., general rate model (GR) as well as pseudo-homogeneous models, i.e., the transport dispersive model (TD) and the equilibrium dispersive model

SC R

(ED). Afterward, we provide a step-by-step guide on how to determine the model

parameters along with their concentration and flow rate dependencies. We exploits the moment analysis method to examine the physical meaning of the model parameters and the manner of their determination, including the effect of protein dispersion in system

U

void volumes and its flowrate dependence. We point out problems that can be

N

encountered when using the moment analysis method. We also describe a methodology

A

for calibration of the model based on several chromatographic runs, which is verified over a wide range of the process conditions pertaining to the column dimensions, mass loading

M

and residence time. As a model system we selected ion-exchange chromatography (IEX), which belongs to one of the most frequently used tools in protein downstream processing.

ED

The analysis is presented for two proteins, i.e., basic fibroblast growth factor 2 (FGF2) and

2. Theory

PT

lysozyme (LYZ), which were eluted from two different IEX media.

2.1. Pseudo-homogenous models: ED and TD

CC E

The ED model is the simplest pseudo-homogenous model which is readily used for

simulating chromatographic elution of small molecules in a chromatographic column of limited efficiency [18]. This model was also used to describe protein chromatography [7,

A

17]. The model consists of only one differential mass balance equation for the liquid phase, which minimizes computational efforts for its solution. The model equation can be expressed as follows: 𝜕𝐶 𝜕𝑡

+𝐹

𝜕𝑞 ∗ 𝜕𝑡

𝑢 𝜕𝐶

+𝜀

𝑡

𝜕𝑥

= 𝐷𝑎,ED

𝜀𝑒 𝜕 2 𝐶 𝜀𝑡 𝜕𝑥 2

(1)

where C [mg mL-1] is the concentration of species in the mobile phase; q* [mg mL−1 bed,solid ] is the concentration in the adsorbed phase at equilibrium (related to the solid matrix 4

volume – mLbed,solid); u [m s-1] is the superficial velocity; t [s], x [m] are the time and axial coordinates, respectively; F is the phase ratio, 𝐹 =

1− 𝑡

𝑡

; t, e and p are the total and

external bed porosities and internal particle porosity, respectively, correlated as follows: 𝑡 = 𝑒 + (1 − 𝑒 ) 𝑝

(2)

The coefficient Da,ED [m2 s-1] in Eq. (1) is the apparent dispersion coefficient of the ED model that lumps the contribution of all non-ideal effects, i.e., axial dispersion and mass

IP T

transport resistances.

The TD model is frequently used for prediction of protein chromatography (e.g., [9, 23,26-29]). It consists of two equations, i.e., the differential mass balance equation in the 𝜕𝐶 𝜕𝑡

+𝐹

𝜕𝑞 𝜕𝑡

𝑢 𝜕𝐶

+𝜀

𝑡

𝜕𝑥

= 𝐷𝑎,TD

SC R

mobile phase that is expressed as follows: 𝜀𝑒 𝜕 2 𝐶 𝜀𝑡 𝜕𝑥 2

(3)

U

and the kinetic equation that describes the linear driving force in the adsorbed or in the 𝜕𝑞

= 𝑘𝑚,𝑞 (𝑞 ∗ − 𝑞) = 𝑘𝑚,𝑐 (𝐶 − 𝐶 ∗ )

(4)

A

𝜕𝑡

N

liquid phase:

M

* -1 where q [mg mL−1 bed,solid ] is the local adsorbed phase concentration; C [mg mL ] is the

equilibrium concentration in the liquid phase; Da,TD [m2 s-1] is the apparent dispersion

ED

−1 coefficient of the TD model, km,q [s-1], and km,c [mL mL−1 bed,solid s ] are the overall lumped

mass transport coefficients that lump external and internal mass transport resistances. The kinetic equations of solid and liquid film driving forces (Eq. (4)) can be used

PT

interchangeably, provided that the coefficients km,q and km,c are properly correlated. When the kinetic parameters km,q and km,c are sufficiently high, the TD model converts to the ED

CC E

model.

The equilibrium concentration in the adsorbed phase, q*, is given by the adsorption

isotherm. The thermodynamics of the IEX process are usually described by the steric mass

A

action model, SMA [30]. When sodium chloride is used as a source of counter ions, which is often the case in IEX chromatography, the SMA equation under non-linear isotherm conditions can be expressed as: 𝑞∗ =

𝐾𝑒 [𝑞∞ − (𝑧+ ) 𝑞 ∗ ]𝑧 [𝑁𝑎+ ]𝑧

𝐶

(5)

where 𝑞 ∞ is the adsorbent capacity; Ke is the equilibrium constant; z is the protein charge;

 is the shielding factor; Na+ is the concentration of counter ions. 5

The implicit form of Eq. (5) with respect to the adsorbed phase concentration can be converted into an explicit function with respect to the liquid phase concentration, as follows: 𝐶∗ =

[𝑁𝑎+ ]𝑧 𝐾𝑒 [𝑞 ∞ − (𝑧+ ) 𝑞 ∗ ]𝑧

𝑞

(6)

For the linear isotherm, the term (z+)q* is neglected, and both Eqs (5) and (6) have an explicit form, i.e.: 𝐻

𝑞

(7)

IP T

1

𝑞 ∗ = 𝐻 𝐶 or 𝐶 ∗ =

where H is the Henry constant expressed as:

SC R

𝐻=

𝐾𝑒 (𝑞 ∞ )𝑧 [𝑁𝑎+ ]𝑧

(8)

The contribution of the kinetic rate of adsorption-desorption to band broadening of

U

elution profiles in IEX chromatography is usually of minor importance compared to the influence brought by mass transport effects; therefore, it was neglected in the analysis.

N

The ED and TD models are combined with the initial conditions defining the initial

A

concentration profile inside the column and the boundary conditions. If the injection

M

profile can be approximated by a rectangular pulse, the boundary conditions at the column inlet (x = 0) and the column outlet (x = L), are expressed as:

𝜕𝐶(𝑡,𝑥 = 𝐿)

=0

(9) (10)

PT

𝜕𝑥

𝐶𝐹 for 𝑡 ∈ [0, 𝑡𝑖𝑛𝑗 ] } 0 for 𝑡 > 𝑡𝑖𝑛𝑗

ED

𝐶𝐹 (𝑡, 𝑥 = 0) = {

2.2. Moment analysis

CC E

2.2.1. Formulation of the HETP equation

The HETP equation originates from the moment analysis method derived based on the

analytical solution of the GR model. For the sake of brevity, the formulation of the GR

A

model is presented in Appendix A. HETP is defined as follows [18]: HETP =

2 𝐿𝑡𝑜𝑡 𝜎𝑡𝑜𝑡

𝜇12

=

𝐿𝑡𝑜𝑡 N

(11)

2 where: 𝜎𝑡𝑜𝑡 is the total peak variance which determines band broadening at the column

outlet; Ltot is the column length; N is the number of theoretical plates; 1 is the first absolute moment, which can be calculated as follows:

6

1 =

∞

∫0 𝑡 𝐶(𝑡) 𝑑𝑡

(12)

∞

∫0 𝐶(𝑡) 𝑑𝑡

If a rectangular injection profile is assumed, then 1 is given by the following equation [18]:

1 =

𝑡𝑖𝑛𝑗 2

𝐿 𝑒

+ 1,𝑣 +

(1 + 𝑘1 )

𝑢

(13)

where the first term of Eq. (13) corresponds to the retention of injection pulse, with tinj

IP T

denoting the injection time; the second term, 1,𝑣 , denotes the peak retention in the total void volume of the chromatographic system, which includes the contribution of the

SC R

external column void volume between the injection system and the detector, and the internal column void volume in the column hardware; the last term corresponds to the peak retention in the packed bed; the coefficient k1 is expressed as follows:

1− 𝜀𝑒 𝜀𝑒

, 𝐹𝑝 =

1− 𝜀𝑝 𝜀𝑝

.

N

where 𝐹𝑒 =

(14)

U

𝑘1 = 𝐹𝑒 𝜀𝑝 (1 + 𝐹𝑝 𝐾)

The coefficient K in Eq. (14) is defined as follows: for the linear isotherm it is equal to

A

the isotherm slope, i.e., the Henry constant (K = H), whereas for the nonlinear isotherm

M

and the breakthrough curve formation it represents the isotherm cord: K = q*/C. For intermediate cases, i.e., chromatographic peaks eluted under non-linear isotherm cord (section 4.2.3.2).

ED

conditions, K can only be approximated by using either the isotherm slope or the isotherm

PT

The peak variance can be expressed as a sum of a few contributions associated with: 2 dispersion in the total system void volume (𝜎𝑎𝑥,𝑣 ), axial dispersion in the packed bed

CC E

2 2 2 (𝜎𝑎𝑥,𝑐 ), shape of the injection profile (𝜎𝑖𝑛𝑗 ), the external mass transport resistances (𝜎𝑒𝑥𝑡 ), 2 and the internal mass transport resistances (𝜎𝑖𝑛𝑡 ), as follows: 2 2 2 2 2 2 𝜎𝑡𝑜𝑡 = 𝜎𝑎𝑥,𝑣 + 𝜎𝑖𝑛𝑗 + 𝜎𝑎𝑥,𝑐 + 𝜎𝑒𝑥𝑡 + 𝜎𝑖𝑛𝑡

(15)

A

2 The value of 𝜎𝑡𝑜𝑡 can be calculated on the basis of the peak second central moment:

2𝑡𝑜𝑡 =

2

∞

∫0 (𝑡 − 1 ) 𝐶(𝑡) 𝑑𝑡

(16)

∞

∫0 𝐶(𝑡) 𝑑𝑡

2 To predict 𝜎𝑡𝑜𝑡 , the moment analysis method can be used. The peak variance

corresponding to the analytical solution of the GR model can be expressed as [18]: 2 2 𝜎𝑡𝑜𝑡,GR = 𝜎𝑎𝑥,𝑣 +

2 𝑡𝑖𝑛𝑗

12

+

2 𝐿 𝐷𝐿,𝑐 𝜀𝑒3 𝑢3

(1 + 𝑘1 )2 +

2 𝐿 𝑘12 𝜀𝑒

𝑑𝑝

𝑢

6 𝐹𝑒

[𝑘

1 𝑒𝑥𝑡

+

𝑑𝑝 10 𝐷𝑒𝑓

]

(17) 7

The first term of Eq. (17) is attributed to dispersion in the total system void volume, the second one – to broadening of a rectangular pulse of the sample, the third one – to axial dispersion in packed bed, the forth one – to the external (kext [m s-1] ) and internal (Def [m2 s-1]) mass transport resistances, where dp [m] is the particle diameter. The coefficient DL,c [m2 s-1], which accounts for axial dispersion in the packed bed, is a linear function of the flow velocity, particle diameter, and molecular diffusion Dm [m2s-1], according to the correlation: (18)

IP T

𝐷𝐿,𝑐 = 𝛾1 𝐷𝑚 + 𝛾2 𝑑𝑝 𝑢

where 1 and 2 are empirical coefficients; Dm is the molecular diffusion coefficient

SC R

(Appendix B, Eq. (B7)).

The coefficients of Eq. (18) are often assessed using empirical correlations (Appendix B, Eqs (B1) and (B2)).

U

The external mass transport coefficient, kext, is a non-linear function of the mobile phase velocity:

N

𝑘𝑒𝑥𝑡 = 𝛾3 𝑢0.333

(19)

A

where 3 is a function of physicochemical properties of the mobile phase, particle

M

diameter, molecular diffusivity, and packing quality, according to literature correlations (Appendix B, Eq. (B3-B7)). The function given by Eq. (19) is concave downward with

ED

respect to u.

The effective diffusion coefficient Def includes the contribution of pore diffusion, and

PT

surface diffusion that can occur in parallel to the pore diffusion: 𝐷𝑒𝑓 = 𝜀𝑝 (𝐷𝑝 + 𝐹𝑝 𝐾 𝐷𝑆,0 )

(20)

CC E

where Dp is the pore diffusion coefficient, which magnitude can be assessed based on empirical correlations (Appendix B, Eqs (B8-B12)); Ds,0 is the surface diffusion coefficient (Appendix B, Eq. (B9)). The flow geometry within void volumes of chromatographic systems is complex and

A

difficult to describe. Therefore, we suggest a simplified way for assessing the protein band 2 broadening in the void volumes (𝜎𝑎𝑥,𝑣 ) by implementing into the model an apparent

column length Lv, which accommodates retention in the total system void volume at the average flow velocity equal to that in the packed column. In turn, the protein dispersion in the total void volume is characterized by an average dispersion coefficient DL,v. Thus, 2 𝜎𝑎𝑥,𝑣 can be expressed as follows:

8

2 𝜎𝑎𝑥,𝑣 =

2 𝐿𝑣 𝐷𝐿,𝑣 𝜀𝑒3

(1 + 𝑘1,0 )

𝑢3

2

(21)

where 𝑘1,0 = 𝜀𝑝 𝐹𝑒 , which corresponds to the coefficient k1 under non-binding conditions 2 (K = 0 in Eq.(14)). Eq. (21) allows evaluating DL,v when 𝜎𝑎𝑥,𝑣 is known.

Consequently, the first absolute moment can be expressed as: 𝑡𝑖𝑛𝑗

1 = 𝐿𝑣 𝑒 𝑢

𝐿 𝑒 𝑢

(1 + 𝑘1 ) +

𝐿𝑣 𝑒 𝑢

(1 + 𝑘1,0 )

(22)

(1 + 𝑘1,0 ) = 1,𝑣 , as the first absolute moment in total void volumes.

IP T

with

2

+

The HETPGR value can be calculated according to Eq. (11), by combining Eqs (17) and (22) at Ltot =L+Lv. 2 𝐿𝑡𝑜𝑡 (𝜎𝑎𝑥,𝑣 +

=

3 𝑡2 2 𝐿 𝑘2 𝑖𝑛𝑗 2 𝐿 𝐷𝐿 𝜀𝑒 1 𝜀𝑒 𝑑𝑝 (1 + 𝑘1 )2 + + 3 12 𝑢 6 𝐹𝑒 𝑢

(

SC R

HETPGR =

𝐿𝑡𝑜𝑡 2𝑡𝑜𝑡 21

𝑡𝑖𝑛𝑗 𝐿  𝐿  + 𝑒 (1+𝑘1 )+ 𝑣 𝑒 (1+𝑘1,0 )) 2 𝑢 𝑢

𝑑𝑝 1 + ]) 𝑘𝑒𝑥𝑡 10 𝐷𝑒𝑓

[

2

(23)

U

In case of non-binding conditions (k1 = k1,0) and a short injection time (i.e., the HETP equation can be derived: 𝐿𝑡𝑜𝑡 ( 𝑒 (1+𝑘1,0 ))

2

+

2 𝐷𝐿,𝑐 𝑒 𝐿 𝑢 𝐿𝑡𝑜𝑡

+2

2 𝑘1,0

𝑢 𝑑𝑝 𝐿

A

2 𝜎𝑎𝑥,𝑣 𝑢2

2 (1+ 𝑘1,0 ) 6 𝜀𝑒 𝐹𝑒 𝐿𝑡𝑜𝑡

[𝑘

1 𝑒𝑥𝑡

+

𝑑𝑝 10 𝐷𝑒𝑓

]

(24)

M

HETPGR =

N

contribution of the variance of the injection profile is negligible), a simplified form of the

Including the flow rate dependencies DL,c and kext, (Eqs 18,19) in Eq. (24), one obtains a

ED

function of HETP versus the mobile phase velocity, u, which represents a modified van Deemter type equation:

B

PT

HETP = HETP𝑣 + A +

𝑢

+ C 𝑢0.667 + D 𝑢

(25)

where HETP𝑣 includes the contribution of dispersion in the void volumes (the first term

CC E

of Eq. (24)); A, B, C, D are the coefficients that can be obtained by matching Eqs (18,19, 24) and (25):

A

A = 2 𝛾2 𝜀𝑒 𝑑𝑝 𝐿

𝐿

B = 2 𝛾1 𝜀𝑒 𝐷𝑚 𝐿 C= 2 D= 2

2 𝑘1,0

(1+𝑘1,0 )

2

2 𝑘1,0

(1+𝑘1,0 )

2

(25a)

𝑡𝑜𝑡

𝐿

(25b)

𝑡𝑜𝑡

𝑑𝑝

𝐿

(25c)

6 𝜀𝑒 𝐹𝑒 𝛾3 𝐿𝑡𝑜𝑡 2 𝑑𝑝

𝐿

60 𝜀𝑒 𝐹𝑒 𝐷𝑒𝑓 𝐿𝑡𝑜𝑡

(25d)

9

When the value of Lv is much smaller than L, HETPv and the term L/Ltot can be neglected and Eq. (25) converts to the common form of the van Deemter equation. Obviously, Eq. (25) is not correct when the injection time is noticeable and its contribution cannot be neglected.

2.2.2. Compatibility conditions for the ED, TD, GR models The moment analysis for the ED and TD models provides adequate moment equations for the ED model: 2 2 𝜎𝑡𝑜𝑡,ED = 𝜎𝑎𝑥,𝑣 +

-

2 𝑡𝑖𝑛𝑗

12

+

2 𝐿 𝐷𝑎,𝐸𝐷 𝜀𝑒3 𝑢3

(1 + 𝑘1 )2

(26)

for the TD model: 2 2 𝜎𝑡𝑜𝑡,TD = 𝜎𝑎𝑥,𝑣 +

2 𝑡𝑖𝑛𝑗

12

+

2 𝐿 𝐷𝑎,𝑇𝐷 𝜀𝑒3 𝑢3

(1 + 𝑘1 )2 +

2 𝐿 𝑘 𝜀𝑡 𝑢

U

where k is defined as:

1 𝑘𝑚,𝑞

SC R

-

IP T

[18, 23-25, 31]:

N

𝑘= 𝐹𝐾

(27)

(28)

A

For the linear isotherm, the coefficient k is equivalent to the retention factor (K = H).

M

The first absolute moment defined by Eq. (13) or (22) is the same for all GR, ED and TD models:

(29)

ED

1,GR = 1,ED = 1,TD

The compatibility conditions for the lumped models is provided by the equality of the

PT

HETP values, i.e.:

(30)

HETPTD = HETPGR

(31)

CC E

HETPED = HETPGR

Since 1 is expressed in the same way for all models, the HETP equality can be

converted to the equality of the second central moments (Eq. (11)):

2𝑡𝑜𝑡,GR = 2𝑡𝑜𝑡,ED = 2𝑡𝑜𝑡,TD

A

(32)

Therefore, the comparison of the left-hand sides of Eqs (17) and (26) allows the determination of the apparent dispersion coefficient for the ED model, Da,ED, provided that the parameters of Eq. (17) are known, i.e.: 𝑘2

𝐷𝑎,ED = 𝐷𝐿,𝑐 + (1+𝑘1

𝑢2 𝑑𝑝

2 6 2 𝐹 1) 𝑒 𝑒

[𝑘

1 𝑒𝑥𝑡

+

𝑑𝑝 10 𝐷𝑒𝑓

]

(33)

10

In the case of the TD model, two coefficients must be determined independently: Da,TD and km. It should be borne in mind that Da,TD is an apparent dispersion coefficient, which is not equivalent to the coefficient DL,c that quantifies axial dispersion in the packed bed. The value of Da,TD can be calculated using the equality given by Eqs (17) and (27) expressed for non-binding conditions, for which K = 0, thus 𝑘1 = 𝑘1,0 = 𝜀𝑝 𝐹𝑒 . Then, one obtains: 2 𝑘1,0

(1+𝑘1,0 )

2

𝑢2 𝑑𝑝 6 𝜀𝑒2 𝐹𝑒

[𝑘

1

+

𝑒𝑥𝑡

𝑑𝑝 10 𝐷𝑒𝑓

]

(34)

IP T

𝐷𝑎,TD = 𝐷𝐿,𝑐 +

The contribution of the second term in Eq. (34) is small due to a low value of 𝑘1,0 , and is

SC R

neglected in literature. However, it plays an important role in case of elution of macromolecules (large mass transport resistances) under non-binding or weakly binding

conditions, as it can determine the column efficiency and its velocity dependence. When (27) for the case of binding conditions: (1+𝑘)2 𝑘

𝜀𝑡 𝑑𝑝 6 𝜀𝑒 𝐹𝑒

𝑘2

((1+ 𝑘1

1

)2

−

2 𝑘1,0

2

(1+ 𝑘1,0 )

) [𝑘

1

𝑒𝑥𝑡

+

𝑑𝑝

10 𝐷𝑒𝑓

]

(35)

A

𝑘𝑚,𝑞

=

N

1

U

Da,TD is quantified, then km can be calculated independently by comparing Eqs (17) and

According to Eq. (35), km is a function of the mobile phase velocity through the kext

M

coefficient (Eq. (19)) and a function of the isotherm shape through k1 (Eq.(14)) and

3. Experimental

ED

through Def (Eq.(20)). The latter occurs when surface diffusion is active in mass transfer.

PT

3.1. Instruments and chemicals 3.1.1. Instruments

Äkta pure 25 (elution of FGF2) and Äkta purifier (elution of LYZ) with UV and

CC E

conductometric detectors and a data station (GE Healthcare Life Sciences, Uppsala, Sweden) were used for chromatographic experiments. The injector was an injection valve with a sample loop.

A

3.1.2. Materials and process conditions

All experiments were performed at room temperature (av. 22°C). Two different

chromatographic systems were used for two model proteins. FGF2 was overexpressed in E. coli and purified via the Heparin Sepharose FF resin (GE Healthcare, method described in Appendix C with FGF purity analysis on SEC). Lysozyme from chicken egg white (LYZ) was purchased from Sigma–Aldrich (Sigma–Aldrich, Poznań, Poland). In the case of the 11

FGF2 system, Tricorn columns were used (GE Healthcare) packed with CM Sepharose FF (GE Healthcare, particle diameter 90 m, pore diameter about 30nm [1]); in the LYZ system, Omnifit columns (Sigma-Aldrich) were used packed with UNOsphere S (Bio-Rad Laboratories, Hercules, CA, USA, particle diameter 80 m, pore diameter about 130 nm [1]). In each system, different process conditions and adsorbents were used. The system characteristics for FGF2 are presented in Table 1.

IP T

The characteristics of the system for LYZ are presented in Table 2. Insert Table 1

SC R

Insert Table 2

3.1.3. Procedures

U

3.1.3.1. Determination of the system void volumes and the bed porosities

To determine the internal void volume (in the column hardware), the empty column of

N

a defined volume and the same column filled with water and closed with stop-plugs were

A

weighed. The difference in the weights was used to calculate the internal void volume for

M

both FGF2 and LYZ systems.

The external void volume for the FGF2 system was determined by injecting pulses of

ED

acetone and protein solutions into the system, in which the column was replaced with a zero-length connector. A flow restrictor was mounted into the system to mimic the presence of the column. The retention time recorded for the “connector” peaks was

PT

calculated based on the first absolute moment and used to determine the external void volume of the system (connection capillaries and valves), which made up a minimum of

CC E

75% (depending on the column size) of the total void volume (external plus internal). To assess the value of the dispersion coefficient in the external void volume, the second central moments of the connector peaks were calculated.

A

In the case of the LYZ system, the external void volume comprised only 18% of the total

void volume, whereas the internal void volume accounted for the remaining 82%. To determine dispersion in total void volume in that system, the pulses of the protein were recorded under non-binding conditions for Col. 1 and Col. 2. The calculation procedure for quantifying dispersion of the protein in the void volumes is described in section 4.3. To determine the total bed porosity in both FGF2 and LYZ systems, the pulses of acetone were injected into the systems equipped with the columns and eluted with a 12

mobile phase of varied salt content: 0-1.5 M NaCl. Additionally, the protein pulses were recorded under non-binding conditions, i.e., at 1.5 M NaCl in the mobile phase. The retention times of the pulses in the column and in void volumes were used to calculate the total bed porosity. The external bed porosity for both systems was determined from the retention time of blue dextran pulses, which were injected into the columns.

IP T

3.1.3.2. Pulse experiments The samples of the protein solutions were injected into the system by the injection

valve and eluted with the mobile phase isocratically, i.e., at a constant salt concentration

SC R

(NaCl) in the solution, or in PB (phosphate buffer free of NaCl). The injection conditions

were: the protein concentration Cp = 1-4 mg mL-1 and the injection volume Vinj = 0.01 mL for low-loading conditions (linear isotherm), and 2-10 mL for high-loading conditions

U

(non-linear isotherm). The protein peaks were recorded at 280 nm wavelength and converted into the concentration units using a calibration factor. For FGF2, the

N

corresponding molar extinction factor of 15930 M-1 cm-1 was calculated by Prot Param

A

program (http://www.expasy.org/tools/protparam.html). In the case of LYZ, the protein

M

standard solutions were used to correlate the detector signal with the protein concentration.

ED

3.1.3.3. Determination of the loading capacity

The loading capacity of CM Sepharose FF resin was determined by the titration method; the value of 𝑞 ∞ = 1.0 [mol L-1 of the bed volume] was obtained. The loading

PT

capacity for UNO-sphere was provided by the manufacturer, and it was equal: 𝑞 ∞ = 1.3 −1 [mol L-1 of the bed volume]. Those values were converted into g L−1 bed,solid or mol Lbed,solid

CC E

(per liter of the solid matrix volume) when used for the model calculations. 3.1.3.4. Determination of the ionic strength of PB The ionic strength of PB buffer was determined by eluting its samples through the

A

conductometric detector of the Äkta system. The detector signal was compared to that obtained for NaCl solutions and the equivalent NaCl concentration for PB buffer was estimated thereof.

13

4. Results and Discussion 4.1. Compatibility of the ED, TD, and GR models for a low-efficiency column A numerical study was performed to evaluate the range of compatibility of the ED, TD, and GR models for low and very low column efficiency. A series of chromatograms were simulated for different column efficiencies and adsorption conditions. The column efficiency was adjusted by manipulating the values of the coefficients: kext, Def , and DL,c.

IP T

The numerical methods used for solving the models were the same as those described in [24] for the GR model, and in [32] for the TD and ED model. The simulations for the ED

and TD models were performed using the lumped parameters Da,ED, Da,TD, km calculated

SC R

according to the compatibility conditions, i.e., Eq. (33) for Da,ED, in the ED model, and Eqs

(34, 35) for the TD model. Typical results of the numerical study are illustrated in Figs 1A and 1B. The analysis of the obtained results revealed that the solutions of the GR and TD

U

models converged even at a very low column efficiency, i.e., for the number of theoretical plates, N, less than 10 (Fig. 1A). Discrepancies between the solutions of both models

N

appeared only for an extremely low column efficiency, N < 5, and at the Biot number Bi >

A

40 (Bi = (kext dp)/Def), at which the contribution of pore diffusion dominated in the mass

M

transport mechanism, while the external mass transport resistances were negligible. In cases of Bi < 40 (data not shown) the solutions of TD and GR models converged even for

ED

N < 5.

Insert Fig. 1

PT

The ED model was found to be much less accurate compared to the TD model. Even though the second central moments of the peaks simulated by both the GR and ED models

CC E

were the same, they differed in the peak skewness related to the third central moment, which was not considered in the compatibility conditions. The differences diminished for N > 20 (Fig.1B).

A

This outcome indicated that the TD model is more suitable for predictions of protein chromatography compared to the ED model, particularly for small scale laboratory columns with a very low efficiency, which are used in high throughput experiments. The TD model provides a good compromise between the model accuracy and the computation time usage required for its solution.

14

4.2. Model building Before the model can be used, the underlying kinetic coefficients must be determined. For this purpose we suggest a procedure consisting of several steps including: column characterization, assessment of protein dispersion in the system void volumes and in the packed bed, the isotherm and kinetic coefficients. The workflow scheme of the procedure is illustrated in Fig. 2.

IP T

Insert Fig. 2

4.3. FGF2 system

SC R

4.3.1. System characterization

The column characterization consisted of the measurements of the packed bed porosities and the system void volumes, including the external (outside the column, in lines between

U

the sample injector and the detector) and internal (in the column hardware) void volumes (section 3.1.3.1). The total system void volume (internal plus external) made up 40% of

N

the bed volume for the smallest reference column, Col. 1, whereas for the largest one, only

A

about 10% (Table 1). The contribution of the internal void volumes was relatively small,

M

i.e., a maximum 4% of the bed volume. The retention pattern of the pulses of the protein and acetone in the external void volumes at different flow rates is illustrated in Fig. 3. A

ED

strong band broadening and deformation of protein peaks can be observed, particularly at lower flow rates. This phenomenon was significantly less pronounced for the acetone pulses. This caused the retention time of the protein and acetone in the void volume to

PT

differ, which stemmed from a difference in the flow behavior of the protein and acetone

CC E

solutions. An explanation of that phenomenon requires a separate analysis. Insert Fig. 3

The total bed porosity was determined from the retention time of a tracer, i.e., acetone,

A

which was eluted from the column with mobile phases of different salt contents (section 3.1.3.1). The bed porosity value obtained did not depend on the salt concentration in the mobile phase, which excluded the possibility of bed shrinking in the presence of salt. Additionally, the retention of protein pulses in the column was measured under nonbinding conditions (section 3.1.3.1). The value of the total bed porosity calculated for the protein pulses was smaller than that for acetone, e.g., for the reference column, Col. 1, the porosity value was equal 0.94 for acetone, and 0.84 for the protein. This difference 15

indicated that a part of the adsorbent pores was inaccessible to the protein. Therefore, for further analysis and calculations, the porosity values determined for the protein pulses were adopted. The external bed porosity was determined for the reference column from the retention time of blue dextran and it equaled e = 0.41, whereas the internal particle porosity was calculated from Eq. (2). The latter value (p = 0.73) was set constant for each column of the FGF2 system, since the same adsorbent was used for the packing of all three

IP T

columns.

4.3.2. Determination of the model parameters under non-binding conditions

Further pulse experiments were performed in the presence of the column under non-

SC R

binding conditions, i.e., with the mobile phase containing 1.5 M sodium chloride, for which FGF2 was not retained. The elution profiles of the pulses were recorded at different mobile phase flow rates.

U

As mentioned in section 2.2.1, for the sake of simplicity, the presence of the void

additional distance: 𝐿𝑣 =

𝑉𝑣𝑜𝑖𝑑 𝑆 𝑡

N

volumes was accounted for indirectly by adequately extending the bed length with an (S is the cross-section area of the column). The distance Lv

A

accommodated the residence time of the solute in the total void volume. The values of

M

2 𝜎𝑎𝑥,𝑣 were calculated by the integration of the void volume peaks, i.e. “connector peaks”

(Eq. 16) recorded at different flow rates and used to determine the average dispersion

ED

coefficients in the external void volume, DL,v (Eq. (21)). The results of the calculations are depicted in Fig. 4. Since the contribution of the internal void volume in the column

PT

hardware to the packed bed volume was small, the effect of band broadening in the column hardware was lumped into the coefficient Da,TD.

CC E

Next, the TD model was used to simulate the column peaks recorded under nonbinding conditions. The column was divided into two parts of different lengths: in the first, the length Lv, the dispersion was determined by the coefficient DL,v, whereas in the second, the length L (corresponding to the bed length), the dispersion was described by

A

the apparent dispersion coefficient in the packed bed, Da,TD. The Henry constant, H, was set equal to zero to mimic non-binding conditions in both parts of the column. To determine Da,TD, the peak fitting method was used. Next, the total first absolute moment of the peak and the total peak variance were calculated based on Eqs (14), (22) and (27) at K = 0, and were used to determine the total HETP value according to Eq. (11). The results of calculations are presented in Fig. 4. 16

Insert Fig. 4 The HETP data obtained were used to evaluate the coefficients of the modified van Deemter curve, i.e., the coefficients: A, B, C, D of Eq. (25), and subsequently to determine underlying model parameters. Since the shape of the curve was distinctly non-linear, the contribution of the term (Cu0.667) in Eq. (25) that accounted for the external mass transport resistance was expected to be important. However, the number of the

IP T

coefficients was too large to estimate them independently based on the shape of the curve.

Even though the term B including the molecular diffusion term in Eq. (25) was neglected,

SC R

which can be justified due to low values of the protein diffusivity, as well as on the basis of evaluations of the Peclet number by the literature correlation (Appendix B, Eq. (B1)), the number of three coefficients was still too high to determine them at an acceptable level of the confidence interval. Nevertheless, the HETP data served for the assessment of the

U

curve intercept to calculate the parameter A (Eq. 25), as well as the dispersion coefficient

N

in the packed bed, DL,c (Eq. 24), and for a first guess of the values of the remaining model

A

parameters, kext, Def. Afterward, to accomplish the model formulation, the pulse

M

experiments under binding condition were used.

4.3.3. Determination of the model parameters under binding conditions

ED

4.3.3.1. Linear isotherm conditions

In three steps, we evaluated the model parameters. In the first step, the pulse experiments were performed under binding conditions corresponding to the linear range

PT

of the isotherm (section 3.1.3.2) at different salt concentrations in the mobile phase ranging from 0 to 1.5 M NaCl, but at the same flow rate (Q = 1 mL min-1). The ionic strength

CC E

of PB buffer, used as the mobile phase, was equivalent to 0.162 M NaCl (section 3.1.3), therefore, for the model calculations that value was added to each molar concentration of NaCl in the solutions.

A

For five pulses that differed markedly in the retention pattern, recorded at: 0.1, 0.14,

0.2, 0.3, 0.5 M NaCl, values of km and H were estimated by peak fitting using the TD model. The value of Da,TD was set according to that determined under non-binding conditions for the corresponding mobile phase velocity. The band profiles along with the simulations are presented in Fig. 5A, the corresponding moment analysis is presented in Table S1A in Supplementary Materials.

17

The H values obtained were correlated with the salt concentration (using Eq. (8)), to determine Ke and z, which yielded Ke = 0.301 (for the NaCl concentration expressed in molar concentration) and z = 5.7. To assess the effect of flow rate on band broadening, pulse experiments were performed by varying the flow rates while maintaining a constant salt concentration, under conditions of strong binding (at Csalt = 0.14 M in the mobile phase). The value of km was estimated for each peak using the TD model. The Henry constant was set as that

IP T

determined in the first step described above. The obtained values were correlated with the mobile phase velocity using Eqs (19) and (35), which contained only two unknown

SC R

parameters: 3, in the dependence kext = 3 u0.333, and Def . The band profiles are depicted in Fig. 5B. the corresponding moment analysis is presented in Table S1B in

N

Insert Fig. 5

U

Supplementary Materials.

A

Once those parameters are evaluated, the compatibility conditions (Eq. (35)) was used again in the subsequent third step, to verify the relationship between the lumped

M

coefficient versus the Henry constant, km = f(H), through the coefficient k1 (Eq. (14)) and Def (Eq. (20). For this purpose, the values of km obtained from peak fitting at different salt

ED

concentrations in the first step were compared to those calculated from the compatibility conditions, Eq. (35). The results were agreed closely under conditions of weak adsorption,

PT

i.e., for Csalt > 0.2 M (H < 10, the retention factor k < 2). At lower salt concentrations, the value of km obtained by peak fitting was higher than that calculated from the compatibility

CC E

conditions. The difference increased with increasing H constant, e.g., for H = 24 the estimated value was already two times higher than that obtained from the compatibility conditions. That difference was attributed to the contribution of surface diffusion, which was enhanced with increasing adsorption strength. The Ds,0 coefficient was adjusted

A

according to the calculated difference by use of Eq. (20). The Ds,0 value was low (Table 3), however, the term containing surface diffusion was active under conditions of very strong adsorption. All the coefficients obtained are summarized in Table 3. Additionally, the contribution of dispersion and mass transport resistances to the total variance value was calculated (Fig. 6). The peak variances in the system void volume were determined by the integration of the connector peaks (Fig. 3). It is evident that at lower 18

values of the mobile phase velocity (below 100 cm h-1) the contribution of dispersion effects, particularly those in the void volumes, dominates over mass transport resistances. It diminishes rapidly as the mobile phase velocity increases, whereas the mass transport resistances become of the major importance. Insert Table 3

IP T

Insert Fig. 6 The value of the axial dispersion coefficient in the packed bed, DL,c, the external mass

transport coefficient, kext, and the pore diffusion coefficient, Def, were also evaluated based

SC R

on the literature correlations given by Eqs (B1-B12, Appendix B). The estimated pore diffusion coefficient and that obtained from the correlations are in relatively good agreement (Table 3), whereas, kext is about four times lower than that obtained from the

U

correlation given by Eqs (B3, B4) (Table 3). It is possible that the latter value incorporated

N

effects other than those arising from external mass transport, and attributed to a specific flow behavior of the protein in the column, which were not included in the model.

A

Therefore, kext should be considered as a lumped model coefficient that incorporates the

M

mobile phase velocity dependence of the column efficiency for the protein, rather than a real internal mass transport coefficient.

ED

4.3.3.2. Non-linear isotherm conditions and the model verification The model with all kinetic parameters was extended to account for non-linear isotherm

PT

conditions. To avoid handling the implicit isotherm equation (Eq. (5)), the liquid driving force equation was used instead of the solid film driving force equation (Eq. (4)). Due to

CC E

dilution of chromatographic peaks, the Henry constant was used to calculate the value K in Eq. (14) and Eq. (20). Then, the lumped coefficients in the kinetic equations (Eq. (4)) can be correlated with each other by the Henry constant: km,c = H km,q. The accuracy of such an approach was confirmed by comparing the solutions obtained by the TD and GR

A

models for different column overloads. Typical results of the comparison are shown in Fig. S1 in Supplementary materials. The model was then used to determine the shielding factor in the isotherm equation, (Eq. (5) or Eq. (6)) under non-linear isotherm conditions. The column load was increased by increasing the injection volume of the sample to obtain a peak with a distinct

19

asymmetry. The shielding factor was found by fitting the model simulations to a selected overloaded peak under conditions of strong adsorption (i.e., Csalt = 0.1 M). The whole model was verified by predicting a few chromatograms obtained for different column overloads and different salt concentrations, which were recorded for the operating conditions different from those used for peak fitting. Typical results of the verification are shown in Fig. S2 (Supplementary Materials). Since the model accounts for the influence of the salt concentration on the protein

IP T

retention, it can be straightforwardly used to predict the course of gradient elution. An

illustration of the prediction accuracy is shown in in Fig. S3 in Supplementary materials.

SC R

A good agreement between experiment and simulations was again achieved, which proved effectiveness of the model for prediction and design of IEX.To build the entire

model, 17 peaks of the proteins were used: 5 – under non-binding conditions, 10 – under binding conditions for the linear isotherm range (5 peaks each to determine the influence

U

of salt and of the flow rate), and 2 for the non-linear range of the isotherm. The mass of

N

the protein standard used for quantifying the model parameters was very small – for 15

A

peaks only 0.01 mL with Cp = 1 or 2 [mg mL-1] of the protein was injected, for non-linear conditions the injection was 2 mL with Cp = 2.5 [mg mL-1]. Four more peaks were used for

M

the model verification.

The model can be easily extended to include co-elution of other components with the Eq. 6 [30].

PT

4.3.4. Scaling up

ED

protein (impurities) by implementing a multicomponent isotherm equation into Eq. 5 or

Finally, we applied the model to predict the process at an increased scale. Two different

CC E

column systems were selected for the analysis: Col. 2, which was the same column as the reference one, but with a threefold increase in length, and Col. 3 with a doubled diameter. The column characteristics are presented in Table 1 in section 3.1.2. The kinetic and

A

thermodynamic parameters were set the same as for the reference column (Col. 1). To account for the difference in the column packing, for both columns Col. 2 and Col. 3, the value of the total bed porosity, t, was adjusted using a single FGF2 peak. The process conditions for recording of that peak were: a strong adsorption, for which the contribution of the void volumes to the peak retention was negligible; a linear range of the isotherm; and the mobile phase velocity within the range of interest. Peak fitting was used to estimate t, which determined the peak retention factor, k, through the phase ratio, 20

F, in the term k = F H. The value of t obtained for Col. 2 differed by 2% compared to the reference column, being t = 0.86 versus t = 0.84 for the latter one. This value corresponded to a 2% error in the determination of the bed porosity, which is below the experimental accuracy; however, it resulted in about 15% error when calculating the phase ratio, F, as well as the retention factor, k. In the case of Col. 3, the t value was 0.85. Additionally, several experimental peaks were recorded for both larger columns to

IP T

verify the whole set of the coefficients under the binding conditions. The profiles were simulated without further adjustments. Typical results of simulations are shown in Fig.7; the corresponding moment analysis is presented in Table S2 in Supplementary materials.

SC R

Small discrepancies between simulations and experimental data visible in Fig. 7 stem from the difference in the packing quality of columns. The predictions reproduce well the first absolute moment (Table S2), and less accurately the peak variances, particularly for

very broad peaks. The latter can be partly explained by errors in baseline correction and

U

peak integration. Nevertheless, the accuracy of predictions of the column efficiency is

N

sufficient.

A

Insert Fig. 7

M

Obviously, for larger columns, for which the ratio of the void volume to the column volume is small, the influence of dispersion in the void volume on the column efficiency

ED

diminishes.

PT

4.4. Lysozyme system

The series of experiments reported above was repeated for LYZ, which was eluted in

CC E

another chromatographic system; the system characteristics is presented in Table 2 in section (3.1.2). The results of experiments and simulations are in principle very similar to those obtained for the FGF2 system, therefore, they are only briefly summarized, with a

A

focus on the differences between them. The main difference between the characteristics of the FGF2 system (Table 1) and the

LYZ system (Table 2) was the contribution of the internal void volumes to the total system void volume. In the LYZ system, the volume of the column hardware made up as much as 82% of the total void volumes. In such a case, the total dispersion in the void volumes could not be determined by replacing the column with a zero-length connector, as it was done when the external void volumes revealed their dominating contribution (FGF2 21

system). For that purpose, the data of pulse experiments were exploited, which were acquired under non-binding conditions for Col. 1 and Col. 2, which referred to the same column, but which was packed with different volumes of adsorbent (Table 2). Then, Da,TD and DL,v were estimated by fitting the model simulations to two peaks recorded for both Col. 1 and Col. 2 at the same flow rates. Further procedures were the same as for the FGF2 system. The results obtained for the reference column were used to generate the van Deemter curve (Fig. 8A), which was used for a first estimate of the model coefficients,

IP T

which were verified on the basis of pulse experiments under binding conditions. A typical comparison of the band profiles vs the salt concentration and the mobile phase velocity

SC R

are depicted in Fig. S4 (Supplementary materials). Insert Fig. 8

The obtained parameters (Table 4) were used to determine the contribution of

U

different non-ideal effects to the total peak variance (Fig. 8B). The calculations shown in

N

Fig 8B indicate that in the system used for the elution of LYZ (Table 2), the dispersion

A

effect in the packed bed was of major importance over a wide range of low mobile phase velocities. This can probably be attributed to non-uniform flow distribution in the LYZ

M

system, which was characterized by a very large volume of the column hardware. The pore diffusion resistances were the most active in the process kinetics for velocities

ED

higher than 120 cm h-1. The remaining conclusions concerning kinetic effects are the same as those drawn for the FGF2 system.

PT

The parameters determined for the reference column were used for scaling up the protein elution in a larger column (Col. 2), and scaling down in a smaller column (Col. 3),

CC E

with only adjustment of the bed porosity. Typical results of the model predictions are illustrated in Fig. S5 (Supplementary Materials). Insert Table 4

A

5. Conclusions

An approach for the design of IEX protein chromatography has been presented, which

consisted of the choice of the dynamic model for the process simulations, and development of the procedure for the determination of the model parameters. The accuracy of three models was compared, i.e., GR, ED, TD models. The simulations using the ED model were found to be inaccurate for chromatographic system efficiency below 22

20 theoretical plates. Hence, it should not be used for the description of protein elution in mini-columns, which are characterized by low or very low plate numbers. The TD model was found to be suitable as a good trade-off between the solution accuracy and the computational effort for its solving. The TD model was then used for simulations of two different IEX chromatographic processes, in which two model proteins were eluted: FGF2 or LYZ. For both cases, the same procedure for the determination of the model parameters was used. The first step of the procedure was the quantification of dispersion in the void

IP T

volumes in both systems. Band broadening in the void volumes was found to be affected by the system hydrodynamics and the type of eluting compound. It was significantly

SC R

pronounced for the protein solutions as compared to that for acetone, and it was different for each of the proteins.

Afterward, the band broadening was determined in the whole system under nonbinding as well as binding conditions, for the linear and non-linear isotherm ranges. The

U

procedure for the model calibration required recording only several peaks with a low

N

material consumption. The approach was successfully verified by changing the scale of

A

the process.

Additionally, the contribution of different kinetic effects to the band broadening was

M

analyzed for both protein systems. In the case of a small column, the contribution of dispersion effects was a factor of major importance, particularly at low mobile phase

ED

velocity and it determined the efficiency of the whole system. Therefore, that contribution has to be distinguished from the other kinetic effects while determining the kinetic

PT

parameters. Moreover, the values of the external mass transport coefficient assessed basing on the elution experiments departed notably from the literature data, which was

CC E

attributed to protein flow behavior. This aspect requires a further analysis. The design procedure developed can be also adopted for different chromatographic

techniques provided that proper adjustment of the isotherm equation is made. To facilitate the procedure for the determination of the model parameters, the statistical

A

methods could be combined with mathematical modelling of the column dynamics.

Acknowledgments Financial support of this work by National Science Center Poland (project DEC2016/22/M/ST8/00193) is gratefully acknowledged. 23

Astrid Dürauer and Alois Jungbauer have received support from the Federal Ministry of Science, Research and Economy (BMWFW), the Federal Ministry of Traffic, Innovation and Technology (BMVIT), the Styrian Business Promotion Agency SFG, the Standortagentur Tirol, the Government of Lower Austria and ZIT – Technology Agency of the City of Vienna through the COMET-Funding Program managed by the Austrian

6. Appendix A. GR Model 𝜕𝐶 𝜕𝑡

𝜕2 𝐶

𝑢 𝜕𝐶

+𝜀

= 𝐷𝐿 𝜕𝑥 2 − 𝐹𝑒 𝑘𝑒𝑥𝑡 (𝐶 − 𝐶(𝑟 = 𝑅𝑝 ))

𝜕𝑥

𝑒

SC R

Differential mass balance equation for all species in the mobile phase:

IP T

Research Promotion Agency FFG.

(A1)

𝜕𝑡

+ 𝐹𝑝

𝜕𝑞 𝜕𝑡

=

𝐷𝑒𝑓 1

𝜕

𝜀𝑝 𝑟 2 𝜕𝑟

𝜕𝐶𝑝

𝑟 2 ( 𝜕𝑟 )

(A2)

N

𝜕 𝐶𝑝

U

Mass balance equation for all species in the stagnant fluid phase in the pores:

A

where Cp is the concentration of species in the stagnant phase in the pores; r is the radial

M

coordinate of the bed particle; the remaining coefficients are defined in section 2.2. The model is coupled with adequate initial and boundary conditions. For Eq. A1 the

ED

boundary conditions are analogous to Eqs (9) and (10), whereas for Eq. A2 boundary conditions at r = 0 (particle center) and r = R (particle radius) are expressed as: 𝜕𝐶(𝑡,𝑟) 𝜕𝑟

= 𝑘𝑒𝑥𝑡,𝑖 (𝐶 − 𝐶𝑝 (𝑡, 𝑟 = 𝑅) )

PT

𝐷𝑒𝑓

𝜕𝐶𝑝 (𝑡,𝑟 = 0)

=0

CC E

𝜕𝑟

(A3) (A4)

B. HETP for the GR model To determine the dispersion coefficient [33] correlation can be used: 𝑢 𝑑𝑝 𝐿 𝜀𝑒

= 0.2 + 0.011 Re0.48

(B1)

A

Pe = 𝐷

where Re =

𝑢  𝑑𝑝 

(B2)

where ρ and η are the density and viscosity of the mobile phase, respectively.

24

The Reynolds number in chromatographic processes is very low, e.g., in this study the maximum value was about 0.13, therefore, the second term in Eq. (B1) could be neglected. This yields a linear flowrate dependence of the dispersion coefficient, in agreement with Eq. (18) (section 2.2.1). The coefficient kext is the external mass transport coefficient, which is a function of flow rate according to the Sherwood number correlations, e.g. [34]: 𝐷𝑚

=

1.09 𝑒

(Re Sc)0.333

(B3)

IP T

𝑘𝑒𝑥𝑡 𝑑𝑝

Sh = or [35]:

Sh = 1.85 (𝐹𝑒 Re Sc)0.333

SC R

(B4)

where Re is given by Eq. (B2) and Sc is defined as: 

Sc =  𝐷

𝑚

(B5)

U

Both Eqs. (B3) and (B4) provide the velocity dependence of kext~u0.333 through the

N

Reynolds number.

Tyn and Gusek [36] suggested the following correlation to calculate Dm : 9.2 10−8 𝑀0.333

(B7)

M

𝑇

=

A

𝐷𝑚 

where Dm is in [cm2 s-1], viscosity  of the protein solution in [Pa s], which should account

ED

for the salt concentration effect, temperature T in [K], molecular weight M in [Da]. The effective diffusion coefficient Def includes the contribution of pore diffusion (Dp) and surface diffusion (Ds):

PT

𝐷𝑒𝑓 = 𝜀𝑝 (𝐷𝑝 + 𝐹𝑝 𝐷𝑆 ) where:

CC E

𝐷𝑆 = 𝐾 𝐷𝑆,0 𝐷𝑝 =

𝐷𝑚 

(B8)

(B9) (B10)

A

where Dm is the molecular diffusion coefficient; Ds,0 is the surface diffusion coefficient; if the driving force of the diffusion is the gradient of the surface concentration, which is often assumed [18], the meaning of K is the same as that used in defining k1 in Eq. (14). The symbol  in Eq. (B10) is the tortuosity factor which can be evaluated using the following correlations [18,19]: 𝜒=

(2− 𝜀𝑝 ) 𝜀𝑝

2

(B11) 25

or 𝜒 = 𝜀𝑝 (1 + 1.5 𝐹𝑝 )

(B12)

In IEX both Dp and Ds,0 are lumped coefficients which include the contribution of the influence of the charge of the solute ions on the mass transport [13].

C. Purification and analysis of FGF2 = 1 cm and L = 11 cm in the following steps:

IP T

Purification of FGF2 was done using the Heparin Sepharose FF resin in a column of ID  loading of the clarified homogenate containing FGF2 (100 - 200 mL), Q = 1 mL min-1, mobile phase – 20mM Tris-HCl with 0.5 M NaCl, pH 7.6 (buffer A);

SC R

 washing: with buffer A (6 CV – column volumes)

 elution by two step gradients of the A and B buffers (buffer B: 20mM Tris-HCl with 2 M NaCl, pH 7):

U

- 0-70% B (5 CV),

N

- 70-100% B (3 CV).

Fractions of the purified FGF2 (collected during the second gradient in the elution step)

A

were pooled and analyzed for the protein concentration (see section 3.1.3.2) and purity.

M

Purity analysis was performed on size exclusion analytical column ACQUITY UPLC BEH125 SEC 1.7 μm, 4.6x150 mm (Waters, USA). The loading was 10 µL of the sample, Q

ED

= 0.3 mL min-1, detection in 214 and 280 nm at 30°C. The FGF2 peak obtained from SEC chromatogram corresponded to av. 97% of the

A

CC E

PT

surface of all peaks present in the chromatogram.

26

References [1] G. Carta, A. Jungbauer, Protein Chromatography: Process Development and Scale-up, WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim, 2010. [2] A.A. Shukla, M.R. Etzel, S. Gadam, Process Scale Bioseparations for the Biopharmaceutical Industry, 2007, CRC Press Taylor & Francis Group. [3] T. Bund, Allelein S., Arunkumar A., Lucey J.A., Etzel M.R., Chromatographic purification

IP T

and characterization of whey protein–dextran glycation products, J. of Chromatogr. A 1244 (2012) 98–105

[4] A. S. Rathore, R. Bhambure, High-Throughput process development: I. Process

SC R

chromatography, in: Protein Downstream Processing, Humana Press, 2014, pp. 29–37.

[5] R. Bhambure, K. Kumar, A. S. Rathore, High-throughput process development for biopharmaceutical drug substances, Trends Biotechnol. 29 (2011) 127–135.

[6] R. Bhambure, A. S. Rathore, Chromatography process development in the quality by

U

design paradigm I: Establishing a high-throughput process development platform as a

A

Biotechnol. Progr. 29 (2013) 403–414.

N

tool for estimating characterization space for anion exchange chromatography step,

M

[7] S. Chhatre, N. J. Titchener-Hooker, Review: microscale methods for high-throughput chromatography development in the pharmaceutical industry, J. Chem. Technol. Biotechnol. 84 (2009) 927–940.

ED

[8] E. J. Close, J. R. Salm, D. G. Bracewell, E. Sorensen, Modelling of industrial 1304–1314.

PT

biopharmaceutical multicomponent chromatography, Chem. Eng. Res. Des. 92 (2014) [9] D. Gétaz, A. Butté, M. Morbidelli, Model-based design space determination of peptide

CC E

chromatographic purification processes, J. Chromatogr. A 1284 (2013) 80–87. [10] W. R. Keller, Steven T. Evans, Gisela Ferreira, David Robbins, Steven M. Cramer, Use of MiniColumns for linear isotherm parameter estimation and prediction of benchtop column performance, J. Chromatogr. A 1418 (2015) 94–102.

A

[11] W. R. Keller, S. T. Evans, G. Ferreira, D. Robbins, S. M. Cramer, Understanding operational system differences for transfer of miniaturized chromatography column data using simulations, J. Chromatogr. A 1515 (2017) 154–163. [12] R. Khalaf, J. Heymann, X. LeSaout, F. Monard, M. Costioli, M. Morbidelli, Model-based high-throughput design of ion exchange protein chromatography, J. Chromatogr. A 1459 (2016) 67–77. 27

[13] A. Wiesel, H. Schmidt-Traub, J. Lenz, J. Strube, Modelling gradient elution of bioactive multicomponent systems in non-linear ion-exchange chromatography, J. Chromatogr. A 1006 (2003) 101–120. [14] B. C.S. To, A. M. Lenhoff, Hydrophobic interaction chromatography of proteins. IV. Protein adsorption capacity and transport in preparative mod, J. Chromatogr. A 1218 (2011) 427–440. [15] A. Osberghaus, S. Hepbildikler, S. Nath, M. Haindl, E. von Lieres, J. Hubbuch, and empiric modeling approaches, J. Chromatogr. A 1237 (2012) 86–95.

IP T

Optimizing a chromatographic three component separation: A comparison of mechanistic

SC R

[16] T. C. Huuk, T. Hahn, A. Osberghaus, J. Hubbuch, Model-based integrated optimization and evaluation of a multi-step ion exchange chromatography, Sep. Purif. Technol. 136 (2014) 207–222.

[17] E. J. Close, J. R. Salm, D. G. Bracewell, E. Sorensen, A model based approach for

N

variability, Chem. Eng. Sci. 116 (2014) 284–295.

U

identifying robust operating conditions for industrial chromatography with process

A

[18] G. Guiochon, A. Felinger, D. Shirazi, A. Katti, Fundamentals of Preparative and Nonlinear Chromatography, Sec. Ed., Academic Press, 2006, pp. 281–345.

M

[19] M. Suzuki, Adsorption Engineering, Elsevier, Amsterdam, 1990. [20] T. Gu, Mathematical Modeling and Scale-up of Liquid Chromatography – with

ED

Application Examples, Sec. Ed., Springer Verlag, Berlin-New York, 2015. [21] K. Westerberg, E. B. Hansen, M. Degerman, T. B. Hansen, B. Nilsson, Model-based

PT

process challenge of an industrial ion-exchange chromatography step, Chem. Eng. Technol. 35 (2012) 183–190.

CC E

[22] J. R. Joseph, A. Sinclair, N. J. Titchener-Hooker, Y. Zhou, A framework for assessing the solutions in chromatographic process design and operation for large-scale manufacture, J. Chem. Technol. Biotechnol. 81 (2006) 1009–1020. [23] K. Kaczmarski, D. Antos, H. Sajonz, P. Sajonz, G. Guiochon, Comparative modeling of

A

breakthrough curves of bovine serum albumin in anion-exchange chromatography, J. Chromatogr A. 925 (2001) 1–17. [24] D. Antos, K. Kaczmarski, A. Seidel-Morgenstern, W. Piątkowski, Concentration dependence of lumped mass transfer coefficients. Linear versus non-linear chromatography and isocratic versus gradient operation, J. Chromatogr. A 1006 (2003) 6176. 28

[25] S. Javeed, S. Qamar, W. Ashraf, G. Warnecke, A. Seidel-Morgenstern, Analysis and numerical investigation of two dynamic models for liquid chromatography, Chem. Eng. Sci. 90 (2013) 17-31 [26] C. K. S. Ng, H. Osuna-Sanchez, E. Valery, E. Sorensen, D. G. Bracewell, Design of high productivity antibody capture by protein A chromatography using an integrated experimental and modeling approach, J. Chromatogr. B 889 (2012) 116-126.

IP T

[27] M. Rüdt, F. Gillet, St. Heege, J. Hitzler, B. Kalbfuss, B. Guélat, Combined Yamamoto

approach for simultaneous estimation of adsorption isotherm and kinetic parameters in

SC R

ion-exchange chromatography, J. Chromatogr. A 1413 (2015) 68–76.

[28] R. Muca, W. Marek, W. Piątkowski, D. Antos, Influence of the sample-solvent on protein retention, mass transfer and unfolding kinetics in hydrophobic interaction chromatography, J. Chromatogr. A 1217 (2010) 2812–2810.

U

[29] R. Muca, W. Marek, S. Woś, W. Piątkowski, D. Antos, Isolation of monoclonal antibody

N

from a Chinese hamster ovary supernatant. II: Dynamics of the integrated separation on

A

ion exchange and hydrophobic interaction chromatography media, J. Chromatogr. A 1305 (2013) 64–75.

M

[30] C. A. Brooks, S. M. Cramer, Steric Mass-Action Ion Exchange: Displacement Profiles and Induced Salt Gradients, AIChE J. 38 (1992) 1969–1978.

ED

[31] L. Lapidus, N. L. Amundsen, Mathematics of Adsorption in Beds. VI. The effect of longitudinal diffusion in ion exchange and chromatographic columns, J. Phys. Chem. 56

PT

(1952) 984–988.

[32] R. Muca, W. Piątkowski, D. Antos, Effects of thermal heterogeneity in hydrophobic

CC E

interaction chromatography, J. Chromatogr. A 1216 (2009) 6716–6727. [33] S. F. Chung, C. Y. Wen, Longitudinal dispersion of liquid flowing through fixed and fluidized beds, AIChE J., 14 (1968) 857–866. [34] E. J. Wilson, C. J. Geankopolis, Liquid Mass Transfer at Very Low Reynolds Numbers

A

in Packed Beds, Ing. Eng. Chem. Fundam., 5 (1966) 9–14. [35] T. Kataoka, H. Yoshida, K Ueyama, Mass transfer in laminar region between liquid and packing material surface in the packed bed, J. Chem.Eng. Jpn., 5 (1972) 132–136. [36] M. T Tyn, T. W. Gusek, Prediction of diffusion coefficients of proteins, Biotechnol. Bioeng., 35 (1990) 327–338.

29

Captions

A

CC E

PT

ED

M

A

N

U

SC R

IP T

Fig. 1. Typical simulations of band profiles by the ED, TD, and GR models. A) the Henry constant: H = 20, and H = 50 at Q = 1 [mL min-1] for N = 4, Bi = 6.4 (kext = 5 x10-6 [m s-1] , Def = 7 x10-11 [m2 s-1], DL = 5.5 x10-6 [m2 s-1]); B) H = 20, Q = 0.5 [mL min-1] for N = 22, Bi = 6 (kext = 2 x10-5 [m s-1] , Def = 3 x10-10 [m2 s-1], DL = 5.5 x10-7 [m2 s-1]). The sample concentration CP,inj = 2 [mg mL-1], Vinj = 0.01 mL. The column characteristics are the same as for the reference column (1) in Table 1.

30

A

CC E

PT

ED

M

A

N

U

SC R

IP T

Fig. 2. The workflow scheme for determining the model parameters

31

A

CC E

PT

ED

M

A

N

U

SC R

IP T

Fig. 3. Chromatograms of FGF2 and acetone pulses in the external void volume (connector peaks) recorded at different flow rates (Q); peak of acetone is shown only for the lowest flow rate Q = 0.1 [mL min-1]

32

M

A

N

U

SC R

IP T

Fig. 4. Dependence of HETP, average dispersion coefficient in the total system void volume, DL,v, and the apparent dispersion coefficient in the packed bed, Da,TD, on the mobile phase velocity for the reference column (Col. 1, Table 1).

A

CC E

PT

ED

Fig. 5. Illustration of band profiles of FGF2. (A) Profiles recorded at different salt contents, Cinj = 1 [mg mL-1], Vinj = 0.01 mL, Q = 1 [mL min-1]; B) profiles recorded at different flow rates: Q = 0.1-2 [mL min-1] (u = 31-612 [cm h-1]), Cinj = 1 [mg mL-1], Vinj = 0.01 mL, Csalt = 0.14 M NaCl, H = 24.

33

M

A

N

U

SC R

IP T

Fig 6. Contribution of the variances corresponding to: dispersion in the packed bed and in the void volumes with i = (ax,c) and i = (ax,v), respectively; internal i = (int) and external i = (ext) mass transport resistances, with respect to the total variance value.

A

CC E

PT

ED

Fig.7. Illustration of the accuracy of the model predictions versus experimental band profiles of FGF2 at different salt content (Csalt = 0.1-1.5 M NaCl) and different column dimensions; A) longer column, Col. 2; B) wider column, Col. 3.

34

A

CC E

PT

ED

M

A

N

U

SC R

IP T

Fig. 8. A) Results of the determination of the column efficiency for the reference column (1) (LYZ – Table 2). A) Values of DL,v in void volumes, the apparent dispersion coefficient, Da,TD, and the course of the van Deemter curve; B) contribution of the peak variances for LYZ corresponding to: dispersion in the packed bed and in the void volumes i = (ax,c) and i = (ax,v), respectively; internal i = (int) and external i = (ext) mass transport resistances, to the total variance value.

35

Table 1. Characteristics of the chromatographic system for FGF2 Model protein – FGF2; Bed – CM Sepharose FF; Mobile phase – 0.1M phosphate buffer (PB); pH 7; 0.03 - 1.5 M NaCl Cinj Vinj Vcol Vvoid ID L Q tot mL

mL

Col. 1a

1.0–4.0

0.01–10

0.94

Col. 2

2.0

0.01–2.0

3.05

Col. 3

2.0

0.01–2.0

4.55

mL = % Vcol 0.37 40% 0.37 12% 0.44 10%

u

cm

cm

[-]

[mL min-1]

[cm h-1]

0.5

4.8

0.84

0.1–2.0

31-612

0.5

15.5

0.86

0.1–1.0

31-306

1.0

5.8

0.85

0.4–4.0

31-306

IP T

[mg mL-1]

A

CC E

PT

ED

M

A

N

U

SC R

where Vinj is the injection volume, Vcol is the column volume, ID and L are the column diameter and length, Q is the mobile phase flow rate. “Col. 1” is the reference column, “Col. 2” is the same column, but packed with a higher volume of the adsorbent, “Col. 3” is another column with a larger diameter. a reference column

36

Table 2. Characteristics of the chromatographic system for LYZ Model protein – LYZ; Bed – UNOsphere S; Mobile phase – 0.1 M phosphate buffer (PB); pH 7; 0.15 –1.0 M NaCl Cinj Vinj Vcol Vvoid ID L Q tot mL

mL

Col. 1a

1.0–2.0

0.01–2.0

3.64

Col. 2

2.0–8.0

0.01–2.0

8.96

Col. 3

2.0–4.0

0.01–2.0

0.88

mL = % Vcol 1.27 35% 1.27 14% 1.15 131%

u

cm

cm

[-]

[mL min-1]

[cm h-1]

1.0

4.63

0.80

0.05–3.0

3.8-229

1.0

11.4

0.82

0.5–3.0

38-229

0.66

2.56

0.76

0.1–2.0

18-351

IP T

[mg mL-1]

A

CC E

PT

ED

M

A

N

U

SC R

“Col. 1” is the reference column, “Col. 2” is the same column, but packed with a higher volume of the adsorbent, “Col. 3” is another column with a smaller diameter. areference column

37

Table 3. Model parameters for the FGF2 system assessed using the TD model (TD) and literature correlations (Appendix B). Thermodynamic parameters

Ke [-]

z



0.301

5.7 Kinetic parameters kext x105 Dp x1011 -1 [m s ] [m2 s-1] 5.32u0.33 TD 3.86 TD 27.0u0.33 b 3.05 d 21.3u0.33 c 5.98 e

6.1

[m2 s-1]

0.68u TD 1.1u a

Ds,0 x1012 [m2 s-1] 1.29 TD

A

CC E

PT

ED

M

A

N

U

SC R

correlation in Appendix B: a(B1, at 1.5 M NaCl); b (B3); c (B4); d (B11); e (B12).

IP T

DL,c x103

38

Table 4. Model parameters for the LYZ system assessed using the TD model (TD) and literature correlations (Appendix B). Thermodynamic parameters

Ke [-]

z



0.014

5.33 Kinetic parameters kext x105 Dp x1011 -1 [m s ] [m2 s-1] 9.01u0.33 TD 2.41 TD 32.4u033 b 2.59 d 0.33 c 25.2u 5.92 e

6.1

[m2 s-1]

1.40u TD 1.00u a

Ds,0 x1012 [m2 s-1] 2.14 TD

A

CC E

PT

ED

M

A

N

U

SC R

correlations in Appendix B: a (B1 at 1 M NaCl); b (B3); c (B4); d (B11); e (B12)

IP T

DL,c x103

39