Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Modeling of controlled-shear affinity filtration using computational fluid dynamics and a novel zonal… Francis, Patrick 2011

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2012_spring_francis_patrick.pdf [ 2.34MB ]
Metadata
JSON: 24-1.0059076.json
JSON-LD: 24-1.0059076-ld.json
RDF/XML (Pretty): 24-1.0059076-rdf.xml
RDF/JSON: 24-1.0059076-rdf.json
Turtle: 24-1.0059076-turtle.txt
N-Triples: 24-1.0059076-rdf-ntriples.txt
Original Record: 24-1.0059076-source.json
Full Text
24-1.0059076-fulltext.txt
Citation
24-1.0059076.ris

Full Text

MODELING OF CONTROLLED-SHEAR AFFINITY FILTRATION USING COMPUTATIONAL FLUID DYNAMICS AND A NOVEL ZONAL RATE MODEL FOR MEMBRANE CHROMATOGRAPHY  by  Patrick Francis  B.A.Sc., The University of British Columbia, 2002  A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in  THE FACULTY OF GRADUATE STUDIES (Chemical and Biological Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  November 2011 © Patrick Francis, 2011  Abstract Controlled-shear affinity filtration (CSAF) is a novel integrated bioprocessing technology that positions a rotor directly above an affinity membrane chromatography column to permit protein capture and purification directly from cell culture. The rotor provides a tunable shear stress at the membrane surface that inhibits membrane fouling and cell cake formation allowing for a uniform filtrate flux that maximizes membrane column performance.  However, the fundamental hydrodynamics and mass transfer kinetics  within the CSAF device are poorly understood and, as a result, the industrial applicability of the technology is limited. A computational fluid dynamic (CFD) model is developed that describes the rotor chamber hydrodynamics of the CSAF device. Once evaluated the model is used to show that a rotor of fixed angle does not provide uniform shear stress at the membrane surface. This results in the need to operate the system at unnecessarily high rotor speeds to reach a required shear stress threshold across the membrane surface, compromising the scale-up of the technology. The CFD model is then used to model design improvements that result in an in silico design of a preparative CSAF device capable of processing industrial feedstocks. To describe mass transfer in stacked-membrane chromatography a novel zonal rate model (ZRM) is presented that improves on existing hold-up volume models. The ZRM radially partitions the membrane stack and external hold-up volumes to better capture nonuniform flow distribution effects. Global fitting of model parameters is first used under non-retention conditions to build and evaluate the appropriate form of the ZRM. Through its careful accounting of transport non-idealities within and external to the membrane stack, the ZRM is then shown to provide, under protein retention conditions, a useful framework for characterizing putative protein binding models, for predicting breakthrough curves and complex elution behavior, and for simulating and scaling separations using membrane chromatography.  ii  By elucidating the intrinsic physical processes ongoing in CSAF the mathematical models presented in this thesis represent essential theoretical tools for the further development of the technology; a technology which has the potential to increase productivity and decrease costs in the downstream processing of biopharmaceuticals.  iii  Preface A version of Chapter 2 of this dissertation has been published as: Francis P, Martinez DM, Taghipour F, Bowen BD, Haynes CA. 2006. Optimizing the rotor design for controlled shear affinity filtration using computational fluid dynamics.  Biotechnol  Bioeng 95:1207-1217. As first author, I was involved in formulating the research design together with the principle investigator C.A. Haynes. In the course of the research I constructed and ran all CFD simulations with the guidance of F. Taghipour; built, with the aid of D.M. Martinez, the necessary program code for the low Reynolds number numerical solution (see Appendix C); performed all PIV experiments; and wrote any additional program code that was required (see Appendix A). In addition, I prepared the manuscript with subsequent editing from all co-authors. A version of Chapter 3 has been published as: Francis P, Haynes CA. 2009. Scale-up of controlled shear affinity filtration using computational fluid dynamics. Biotechnol J 4: 665-673. As first author, I formulated the research design together with the principle investigator C.A. Haynes. I also performed all model simulations and prepared the manuscript with subsequent editing from C.A. Haynes. A version of Chapter 4 has been published as: Francis P, von Lieres E, Haynes CA. 2011. Zonal rate model for stacked membrane chromatography. I: characterizing solute dispersion under flow-through conditions. J Chromatogr A 1218:5071-5078. Together with E. von Lieres and C.A. Haynes I instigated and led development of the theoretical framework of the zonal rate model. E. von Lieres assisted in developing the solution strategy and computational algorithm, and construction of the necessary program code (see Appendix E). I generated all of the simulation results described in this chapter as well as performed all experiments. In addition, I constructed and ran supplementary program codes that were required above and beyond the zonal rate model itself (see Appendix F). Finally, I prepared the manuscript with subsequent editing from E. von Lieres and C.A. Haynes.  iv  A version of Chapter 5 has been submitted for publication as: Francis P, Haynes CA, von Lieres E. Zonal rate model for stacked membrane chromatography. II: characterizing ion-exchange membrane chromatography under protein retention conditions.  As  mentioned above E. von Lieres assisted me in constructing the zonal rate model. I implemented the spreading and bi-Langmuir isotherms that were used in this chapter. I also constructed and ran all supplementary program code that was required (see Appendix G). As first author, I prepared the manuscript with subsequent editing from E. von Lieres and C.A. Haynes.  v  Table of Contents Abstract…. ......................................................................................................................... ii Preface…… ....................................................................................................................... iv Table of Contents ............................................................................................................. vi List of Tables .................................................................................................................... ix List of Figures ................................................................................................................... xi Nomenclature ................................................................................................................ xvii Acknowledgements ...................................................................................................... xxiii Chapter 1 Introduction ................................................................................................... 1 1.1 Thesis overview ....................................................................................... 1 1.2 Current techniques in downstream processing of biopharmaceuticals .... 5 1.2.1 Cell separation techniques .............................................................. 5 1.2.1.1 Sedimentation .................................................................... 6 1.2.1.2 Centrifugation .................................................................... 6 1.2.1.3 Dynamic microfiltration .................................................... 7 1.2.2 Protein purification techniques ....................................................... 8 1.2.2.1 Classic packed-bed chromatography ................................. 9 1.2.2.2 Membrane chromatography ............................................ 12 1.2.3 Integrated cell separation and product capture techniques ........... 14 1.2.3.1 Aqueous two-phase extraction ........................................ 14 1.2.3.2 Expanded bed adsorption ................................................ 15 1.2.3.3 High-gradient magnetic separation ................................. 17 1.2.3.4 Controlled-shear affinity filtration .................................. 17 1.3 Thesis objectives .................................................................................... 20 1.4 Computational fluid dynamic modeling ................................................ 21 1.5 Modeling of membrane chromatography............................................... 25 1.6 Thesis outline ......................................................................................... 30  vi  Chapter 2 Optimizing the rotor design for controlled-shear affinity filtration using computational fluid dynamics .................................................................... 32 2.1 Synopsis ................................................................................................. 32 2.2 Computational and experimental materials and methods ...................... 33 2.2.1 Computational fluid dynamics ...................................................... 33 2.2.2 Particle image velocimetry ........................................................... 35 2.3 Numerical solution for low Reynolds numbers ..................................... 36 2.4 Results and discussions .......................................................................... 40 2.4.1 Evaluation of CFD simulations..................................................... 40 2.4.2 CFD modeling of original CSAF device ...................................... 44 2.4.3 Optimization of CSAF rotor geometry ......................................... 46 2.5 Concluding remarks ............................................................................... 50 Chapter 3 Scale-up of controlled-shear affinity filtration using computational fluid dynamics ....................................................................................................... 51 3.1 Synopsis ................................................................................................. 51 3.2 CFD model and computational methods................................................ 52 3.3 Results and discussions .......................................................................... 53 3.4 Concluding remarks ............................................................................... 64 Chapter 4 Characterizing solute dispersion in stacked membrane chromatography under flow-through conditions using a zonal rate model ........................ 65 4.1 Synopsis ................................................................................................. 65 4.2 Theory… ................................................................................................ 66 4.2.1 ZRM model description ................................................................ 66 4.2.2 Computational methods ................................................................ 70 4.3 Experimental materials and methods ..................................................... 72 4.4 Results and discussions .......................................................................... 72 4.4.1 Comparison of existing hold-up volume models with experiment 72 4.4.2 ZRM model configuration optimization and application.............. 76 4.4.3 Decomposition of system dispersion using the ZRM ................... 79 4.5 Concluding remarks ............................................................................... 84  vii  Chapter 5 Characterizing stacked membrane chromatography under protein retention conditions using a zonal rate model .......................................... 86 5.1 Synopsis ................................................................................................. 86 5.2 Theory… ................................................................................................ 87 5.2.1 Zonal rate model .......................................................................... 87 5.2.2 Intrinsic binding rate models ....................................................... 89 5.2.3 Computational methods ............................................................... 92 5.3 Experimental materials and methods ..................................................... 93 5.4 Results and discussions .......................................................................... 94 5.4.1 Dimensionless group analysis ...................................................... 94 5.4.2 Benchmark study – the Thomas model ........................................ 96 5.4.3 Accounting for external flow non-idealities using the ZRM ..... 100 5.4.4 Evaluating putative binding models........................................... 102 5.4.4.1 Langmuir model ............................................................ 102 5.4.4.2 Bi-Langmuir and SMA models ..................................... 103 5.4.4.3 Spreading model ............................................................ 106 5.2 Concluding remarks ............................................................................. 111 Chapter 6 Conclusion .................................................................................................. 113 References. ..................................................................................................................... 119 Appendices ..................................................................................................................... 136 Appendix A:  Description and program code for the determination of shear profiles from Fluent output ................................................ 136  Appendix B:  Description of the similarity solution for a two-domain model in the absence of suction through the bottom wall ............ 141  Appendix C:  Description and program code for the one-domain similarity solution............................................................................... 147  Appendix D:  Description of modeling equations for a CSTR with multiple inflows ............................................................................... 153  Appendix E:  Program code for the zonal rate model .............................. 155  Appendix F:  Program code for the moment analysis of breakthrough curves ................................................................................. 183  Appendix G:  Program code for the determination of Thomas model parameters .......................................................................... 186  viii  List of Tables Table 2.1  Characteristics of the CSAF device used in the CFD model. ................... 34  Table 3.1  Characteristics of the scaled CSAF devices used in the CFD model. ...... 53  Table 4.1  Sum of squared residuals for the fitting of various models to experimental breakthrough data at four different flow rates for loading of ovalbumin (1 g L-1) on a Mustang XT5 column under non-binding conditions. Dx was 7 x 10-11 m2 s-1 for all model calculations. ................................................... 75  Table 4.2  Parameters for the symmetric two-zone ZRM configuration regressed from experimental breakthrough data at four different flow rates for ovalbumin (co = 1 g L-1) loaded on a Mustang XT5 column under nonbinding conditions. Dx was set to 7 x 10-11 m2 s-1. ................................... 78  Table 4.3  Breakdown of ZRM predicted solute elution variance v (s) for the 1.5 mL min-1 breakthrough data as a function of flow zone and flow elements within the zone. ......................................................................................... 81  Table 5.1  Structural parameters for (i) the symmetric one-zone model proposed by Roper and Lightfoot [136] comprised of a PFR in series with the membrane stack sandwiched between two CSTRs of equal residence time and (ii) the symmetric two-zone nr-ZRM. Model parameters were determined by fitting experimental data for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 column under non-binding conditions at flow rates of 1.5, 5, 10 and 20 mL min-1. In both models, all dispersion occurring in the dead volume modeled by the PFR is accounted for by increasing the solute residence time of the first CSTR in the network by a flow rate independent amount ′, so that 1/1 = 1/′ + Q1/VCSTR..................................................................................................... 99  Table 5.2  Regressed parameter values for the Langmuir adsorption model determined from fitting BTC data for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 column under retained conditions at four different flow rates (1.5, 5, 10 and 20 mL min-1). The Langmuir adsorption model was used within the symmetric two-zone ZRM. The rate constants in parentheses represent the Langmuir kinetics parameters determined from fitting the corresponding Thomas model to the BTC data at 20 mL min-1. ......................................................................................................... 99  Table 5.3  Sum of squared residuals for the fitting of various models to BTC data. Experimental conditions are the same as in Table 5.2. The performance of each binding model is that when implemented within the symmetric twozone ZRM. .............................................................................................. 104 ix  Table 5.4  Regressed parameters for the bi-Langmuir adsorption model when introduced into the two-zone ZRM. Experimental conditions and regression method are the same as in Table 5.2...................................... 105  Table 5.5  Regressed parameters for the simplified spreading model when introduced into the two-zone ZRM. Experimental conditions and regression method are the same as in Table 5.2. ................................................................... 107  x  List of Figures Figure 1.1  Comparison of a conventional protein purification stream (left side) to a corresponding bioprocessing scheme where the primary isolation and product capture stages have been integrated (right side). ........................... 3  Figure 1.2  Permeate flux characteristics of conventional tangential flow filtration. The transmembrane pressure decreases with position resulting in an inhomogeneous permeate flux. ................................................................... 8  Figure 1.3  Separation of three solute molecules in elution chromatography. As the mobile phase flows through the column from left to right the three molecules (represented as squares, circles and diamonds) are retained differentially by the stationary phase and thus exit the column at different times. ......................................................................................................... 10  Figure 1.4  Solute transport processes in packed bed chromatography (left side) and membrane chromatography (right side). In packed bed chromatography solute molecules must, in order to reach the adsorption surface on the interior of porous beads, diffuse through a thin film layer surrounding the bead and then the stagnant fluid within the pores. Conversely, in membrane chromatography the solute molecules are convected past the binding sites and the only diffusive process is across the thin film layer surrounding the pore wall. ........................................................................ 12  Figure 1.5  Comparison of a conventional packed bed chromatography column (left side) to expanded bed adsorption (right side). In conventional packed bed chromatography the bed is constrained and would become plugged with the introduction of particulate matter. In expanded bed adsorption the upward flow of mobile phase through an unconstrained bed of adsorbent beads causes the bed to expand creating spaces between the beads that allow for the passage of particulate matter. .............................................. 16  Figure 1.6  Schematic of CSAF. Feed is introduced to the rotor chamber via A and retentate is removed via B. A conically shaped rotor, E, spins above an affinity membrane stack, D. The shear caused by the rotor inhibits cell deposition allowing for uniform filtrate flux, C........................................ 18  Figure 1.7  Schematic of a connected CSAF configuration in which the rotor, E, spins above a microfiltration membrane. The filtrate from this membrane is then passed directly to an external affinity membrane stack, D. Feed is introduced to the rotor chamber via A, retentate is removed via B and filtrate exits the membrane stack via C. .................................................... 19  xi  Figure 1.8  Diagram of a system dispersion model for stacked membrane chromatography wherein the plug flow reactor (PFR) accounts for a lag time uncoupled from system dispersion and the constant stirred tank reactor (CSTR) describes dispersion effects in the extra column hold up volumes (e.g. tubing, valves, etc.). ........................................................... 28  Figure 1.9  Representative path lengths for solute flow within a stacked-membrane chromatography module possessing both axial and radial flow components within the feed distribution and eluent collection manifolds. Each shaded area indicates a radially defined zone in the manifold offering a distinct solute residence time, with the length of the flow line through the zone reflecting the time. Qi represents the flow rate through zone i. ............... 29  Figure 2.1  The tangential component of velocity profiles 1 and 2 estimated in the limit of = 0. 1 and 2 represent the dimensionless tangential velocity components in the porous and fluid domains respectively. Four curves are shown: (a) Da/m = 1 x 10-5, (b) Da/m = 1 x 10-3, (c) Da/m = 1 x 10-2 and (d) Da/m = 1 x 105. These estimates were made with B = 0.15,  = 0.25, H = 2 x 10-3 m and m = 0.58. ....................................................................... 37  Figure 2.2  Physical configuration of the 1-domain low Reynolds number numerical solution of the CSAF hydrodynamics. A conical rotor with angle of incidence  and a localized height above the membrane of h(r) is rotated at a frequency . Fluid flows out through the membrane at a constant velocity Um. ............................................................................................... 38  Figure 2.3  CFD computed () dimensionless (A) tangential (), (B) radial (g’) and (C) axial (g) velocity profiles under the rotor at a radial distance from the rotor center-point of  = 0.71, compared to those calculated by the 1domain similarity solution (solid line). Rrot = 35 mm,  = 10 rpm and Um = 3.65 x 10-5 m s-1. .................................................................................... 41  Figure 2.4  CFD computed () shear stress profile along the membrane surface compared to that calculated by the 1-domain similarity solution (solid line) when Rrot = 35 mm,  = 10 rpm and Um = 3.65 x 10-5 m s-1. ................... 42  Figure 2.5  Tangential velocity profiles determined by PIV measurement () compared to CFD results (▲). Experiments/simulations were carried out with an angled rotor at a rotational frequency of 700 rpm for radial positions of (A)  = 0.57 and (B)  = 0.86 when  = 0. The error bars represent the standard deviation................................................................ 43  Figure 2.6  Shear stress profiles along the membrane surface for the original CSAF rotor at  = 700 rpm () and for the CFD optimized rotor at  = 250 rpm (▲) when Rrot = 35 mm and filtrate flux is 125 L m-2 hr-1. ...................... 45  xii  Figure 2.7  Tangential velocity profiles under the original CSAF rotor at various radial positions: (a)  = 0.29, (b)  = 0.57, (c)  = 0.71 and (d)  = 0.86 when  = 700 rpm, Rrot = 35 mm and the filtrate flux is 125 L m-2 hr-1. . 45  Figure 2.8  Schematic of (a) the CFD optimized rotor shape compared to (b) the original CSAF rotor having a fixed angle of 4o. ....................................... 47  Figure 2.9  Tangential velocity profiles under the CFD optimized rotor at various radial positions: (a)  = 0.29, (b)  = 0.57, (c)  = 0.71 and (d)  = 0.86 when  = 250 rpm, Rrot = 35 mm and the filtrate flux is 125 L m-2 hr-1. . 48  Figure 2.10  Shear stress profiles along the membrane surface for the CFD optimized rotor at different rotor speeds: 200 rpm (), 300 rpm (▲), 400 rpm (), 500 rpm (+) and 700 rpm (). Rrot = 35 mm and filtrate flux is 125 L m -2 hr-1 . ................................................................................................. 49  Figure 3.1  Schematic of a scale-up appropriate CSAF device. A and B refer to the feed and retentate lines respectively and have been shifted from the top of the unit to the sides. C refers to the filtrate from the membrane stack, D. The rotor, E, is supported by a shaft that now extends through the entire CSAF device. ............................................................................................ 52  Figure 3.2  Shear stress profile along the membrane surface for the optimized, smallscale CSAF device.  = 250 rpm, Rrot = 3.5 cm...................................... 54  Figure 3.3  Shear stress profile along the membrane surface for a small-scale CSAF device that has side feed and retentate inlets and a shaft extending through the entire device.  = 250 rpm, Rrot = 3.5 cm and Rshaft = 0.25 cm. ........ 55  Figure 3.4  Shear stress profile along the membrane surface for a CSAF device that has been scaled-up by a factor of 10 fold.  = 250 rpm, Rrot = 35 cm and Rshaft = 2.5 cm. The inset shows that the minimum shear threshold of 0.17 Pa (dotted line) is being maintained over an inner domain that stretches from the shaft to a radial position of  = 0.15. .......................................... 57  Figure 3.5  Dimensionless tangential velocity profiles under a scaled-up CSAF rotor at various radial positions: (a)  = 0.11, (b)  = 0.43 and (c)  = 0.80 when  = 250 rpm and Rrot = 35 cm. ................................................................. 57  Figure 3.6  Shear stress profile along the membrane surface for a CSAF device that has been scaled-up by a factor of 40 fold.  = 20 rpm, Rrot = 140 cm and Rshaft = 27 cm. ............................................................................................ 60  Figure 3.7  Normalized absolute pressure profile along the membrane surface for a CSAF device that has been scaled-up by a factor of 40 fold.  = 20 rpm, Rrot = 140 cm and Rshaft = 27 cm.  = P(r)/Pmax where P(r) is the absolute  xiii  pressure and Pmax is the maximum absolute pressure seen at the membrane surface. ...................................................................................................... 61 Figure 3.8  Normalized axial velocity profile in the affinity membrane stack for a CSAF device that has been scaled-up by a factor of 40 fold.  = 20 rpm, Rrot = 140 cm and Rshaft = 27 cm. ax = uax(r)/uax,max where uax(r) is the axial velocity in the membrane stack and uax max is the maximum axial velocity in the affinity membrane stack. ................................................... 62  Figure 3.9  Proposed configuration for the parallelization of CSAF units. Individual CSAF units are stacked one atop the other with one centralized shaft and individual feed (A), retentate (B) and filtrate (C) lines. ........................... 63  Figure 4.1  Stirred tank and membrane network representation of a stacked-membrane chromatography module with (A) standard and (B) diagonal flow distribution and collection......................................................................... 67  Figure 4.2  Model fit of experimental breakthrough curves (●) for loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 module under non-binding conditions at a flow rate of 1.5 mL min-1.................................................. 74  Figure 4.3  Symmetric two-zone ZRM fit (solid line) of the experimental breakthrough data (●) for ovalbumin (co = 1 g L-1) loaded on a Mustang XT5 module under non-binding conditions at a flow rate of (A) 1.5 mL min-1 and (B) 5 mL min-1. Dx was 7 x 10-11 m2 s-1 in all calculations. ..... 77  Figure 4.4  Contributions of the inner and outer radial zones to elution band spreading of ovalbumin (co = 1 g L-1) loaded on the XT5 module at (A) 1.5 mL min-1 or (B) 20 mL min-1. The symmetric two-zone configuration of the ZRM was used for all calculations. The residence times of the inner and outer tanks were 45.4 and 125.3 s, respectively, and the time lag was 15.25 s. Dx was set to 7 x 10-11 m2 s-1. .................................................................... 80  Figure 4.5  Normalized second moment of ZRM simulated breakthrough curves as a function of axial Peclet number Peax and flow fraction through the outer zone (2) for the case where the column length is 2.2 x 10-3 m and the flow rate is (A) 1.5 mL min-1 or (B) 20 mL min-1. The symmetric twozone configuration of the ZRM was used for all calculations. ................. 82  Figure 4.6  Second central moment of ZRM simulated breakthrough curves as a function of column length and flow fraction through the outer zone (2) for the case where the flow rate is 1.5 mL min -1 and Dx is 7 x 10 -11 m2 s-1…. .................................................................................................... 84  Figure 5.1  Stirred tank and membrane network representation of a stacked-membrane chromatography module with two radial zones. The total flow rate, Q, is split between zone 1, Q1, and zone 2, Q2. The stirred tank network is  xiv  symmetric in that the residence times, i, of the two tanks in each radial zone are equivalent. .................................................................................. 88 Figure 5.2  Schematic of the mechanism of adsorption of a protein molecule according to the simplified spreading model, where protein molecules can reversibly adsorb to a homogeneous binding surface in two different conformations/orientations. ...................................................................... 91  Figure 5.3  Calculated height equivalent to a theoretical plate (HETP), estimated film Peclet number (Pef) and calculated axial Peclet number (Peax) for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 column under binding conditions. HETP values are based on plate theory and derived from breakthrough curves at 1.5, 5, 10 and 20 mL min-1. Film diffusion coefficients are estimated from Equation 5-16 .......................... 96  Figure 5.4  Comparison of the analytical solution of the Thomas model coupled with a pre- and post-column CSTR to experimental breakthrough data (♦) for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 module under binding conditions at flow rates of (A) 1.5 mL min -1and (B) 20 mL min-1. .................................................................................................. 98  Figure 5.5  Model fits of experimental breakthrough data (♦) for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 module under non-binding conditions at flow rates of (A) 1.5 mL min-1 and (B) 20 mL min-1 ....... 101  Figure 5.6  Comparison of the symmetric two-zone ZRM incorporating the Langmuir kinetic model to BTC data for ovalbumin (co = 1 g L-1) loaded on a Mustang XT5 module under binding conditions at flow rates of: 1.5 mL min-1 (■), 5 mL min-1 (▲), 10 mL min-1 (♦) and 20 mL min-1 (●) ......... 103  Figure 5.7  Comparison of the symmetric two-zone ZRM incorporating the biLangmuir kinetic model to BTC data for ovalbumin (co = 1 g L-1) loaded on a Mustang XT5 module under binding conditions at flow rates of: 1.5 mL min-1 (■), 5 mL min-1 (▲), 10 mL min-1 (♦) and 20 mL min-1 (●). . 105  Figure 5.8  Comparison of the symmetric two-zone ZRM incorporating the simplified spreading model to BTC data for ovalbumin (co = 1 g L-1) loaded on a Mustang XT5 module under binding conditions at flow rates of: 1.5 mL min-1 (■), 5 mL min-1 (▲), 10 mL min-1 (♦) and 20 mL min-1 (●) ......... 107  Figure 5.9  Difference in the stability of bound protein states normalized by the projected surface area per bound molecule (G1 – G2/) for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 column at flow rates of 1.5, 5, 10 and 20 mL min-1. ........................................................ 109  Figure 5.10  Experimentally derived dynamic binding capacities, q*, (■) for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 column at flow rates of 1.5, 5, 10 and 20 mL min-1. Two-zone ZRM/spreading model xv  predictions of q*, the concentrations of bound protein in state 1 (q1) and state 2 (q2), and the theoretical dynamic capacity based on all protein binding in state 1 ( q m ,1 ) are also shown. The dashed line represents the saturation capacity (qm) of the membrane for ovalbumin determined from (near-)static isotherm experiments. ......................................................... 110 Figure B.1  Physical configuration of the 2-domain similarity solution for CSAF hydrodyamics. A conical rotor with an angle of incidence  and a localized height, h(r), above the porous medium of thickness Ho, is rotated at a frequency . The axial velocity at z = 0 is specified as Um. ........... 141  Figure D.1  Junction of flows through three zones in a cascade of CSTRs……..154  xvi  Nomenclature av  surface area to volume ratio  [m-1]  An  membrane cross sectional area  [m2]  b  scaling factor  [-]  B  shear stress proportionality constant  [-]  c  solute concentration in mobile phase  [g L-1]  cA  salt ion concentration in mobile phase  [g L-1]  cf  solute concentration in film  [g L-1]  cin  solute concentration at column inlet  [g L-1]  co  initial solute concentration  [g L-1]  cout  solute concentration in column outlet  [g L-1]  cCSTRi  solute concentration in CSTR i  [g L-1]  cinCSTRi  solute concentration in inlet to CSTR i  [g L-1]  c inPFR  solute concentration in inlet to PFR  [g L-1]  solute concentration in outlet from PFR  [g L-1]  PFR cout c~  dimensionless solute concentration in mobile phase  [-]  dp  pore diameter  [m]  Dm  molecular diffusivity  [m2 s-1]  Dx  axial dispersion coefficient  [m2 s-1]  Da  Darcy number  [-]  Dkf  film Damköhler number  [-]  gi  assumed axial velocity function in domain i  [-]  G  change in Gibbs free energy  h  height between rotor and porous domain  [m]  hmax  maximum rotor height from porous domain  [m]  H  total height of rotor from bottom of porous domain  [m]  Hcham  rotor chamber height  [m]  [J mol-1]  xvii  Ho  thickness of the membrane stack  [m]  HETP  height equivalent to a theoretical plate  [m]  k  turbulent kinetic energy  [J] -1 -1  k12  rate constant for spreading from state 1 to state 2  k21  rate constant for unspreading from state 2 to state 1  ka  adsorption rate constant  [L g-1 s-1]  ka,i  adsorption rate constant for site/state i  [L g-1 s-1]  kd  desorption rate constant  [s-1]  kd,i  desorption rate constant for site/state i  [s-1]  kf  film mass transfer coefficient  k f  lumped film mass transfer coefficient  [s-1]  K  membrane permeability  [m2]  K1a  adsorption equilibrium constant  [M-1]  K12  spreading equilibrium constant  [M -1]  L  column length  n  number of zones  P  pressure  [Pa]  Pmax  maximum absolute pressure at membrane surface  [Pa]  Peax  axial Peclet number  [-]  Pef  film Peclet number  [-]  q  solute-ligand concentration in stationary phase  [g L-1]  qi  solute-ligand i concentration in stationary phase  [g L-1]  qm  maximum solute-ligand capacity in stationary phase  [g L-1]  qm,i q~  maximum solute-ligand i capacity in stationary phase  [g L-1]  qA  bound salt co-ion concentration available for exchange  Q  overall volumetric flow rate  [mL min-1]  Qi  volumetric flow rate in zone i  [mL min-1]  r  radial position  [m]  rp  pore radius  [m]  dimensionless solute-ligand concentration in stationary phase  [L g s ] [s-1]  [m s-1]  [m] [-]  [-] [g L-1]  xviii  [J K-1 mol-1]  R  gas constant  Rcham  rotor chamber radius  [m]  Rin  inlet line radius  [m]  Rmem  membrane stack radius  [m]  Rpos  radial position of inlet and retentate lines  [m]  Rret  retentate line radius  [m]  Rrot  rotor radius  [m]  Rshaft  shaft radius  [m]  Re  Reynolds number  [-]  Rey  turbulent Reynolds number  [-]  s  solute concentration in stationary phase  t  time  [s]  tconv  characteristic time for convective mass transfer  [s]  tfilm  characteristic time for film mass transfer  [s]  tlag  lag time  [s]  tplate  characteristic time for convective mass transfer through a single theoretical mass-transfer unit  [s]  ~t  dimensionless time with respect to HETP  [-]  tˆ  dimensionless time with respect to film mass transfer  [-]  T  temperature  u  fluid velocity in porous domain  [m s-1]  uax  axial fluid velocity in porous domain  [m s-1]  uax,max  maximum axial fluid velocity in porous domain  [m s-1]  uint  interstitial fluid velocity  [m s-1]  usup  superficial fluid velocity  [m s-1]  Um  flux through membrane bottom  [m s-1]  v  fluid velocity in rotor chamber  [m s-1]  vr  radial fluid velocity in rotor chamber  [m s-1]  vrot  rotor velocity  [m s-1]  v  tangential fluid velocity in rotor chamber  [m s-1]  v  characteristic charge  [g L-1]  [K]  [-]  xix  Vcolumn  membrane stack volume  [mL]  VCSTR  CSTR volume  [mL]  VPFR  PFR volume  [mL]  +  y  dimensionless wall distance  z ~z  axial position  *  signifies that adsorbed quantity is dynamic  dimensionless axial position  Greek Symbols  [-] [m] [-]      Reynolds number of flow through membrane  [-]    rotational Reynolds number  [-]    rotor gap height    ratio of projected surface areas of spread and unspread protein    turbulent kinetic energy dissipation  f  fraction of stationary phase volume occupied by film  [-]  m  membrane stack porosity  [-]    non-dimensionalized interface height  [-]    non-dimensionalized axial position  [-]    tangential position  i  non-dimensionalized tangential velocity function in domain i    ion capacity in stationary phase    viscosity  [kg m-1 s-1]  t  turbulent viscosity  [kg m-1 s-1]  0  zeroth moment  1   first absolute moment  [s]  2  second central moment  [s2]    kinematic viscosity  ax  normalized axial velocity  [-]    non-dimensionalized radial position  [-]  [m] [-] [J s-1 kg-1]  [rad] [-] [g L-1]  [s g L-1]  [m2 s-1]  xx    normalized absolute pressure    density    steric factor  [-]  v    variance  [s]  [-] [kg m-3]  shear stress  [Pa]   w  shear stress at pore wall  [Pa]    solute residence time in CSTR  [s]  ′  flow rate independent contribution to dispersion  [s]    angle of incidence between rotor and membrane  [rad]  i  flow fraction through zone i  i  streamline function in domain i    rotational frequency  Abbreviations ATP  aqueous two-phase  BHK  baby hamster kidney  BTC  breakthrough curve  CFD  computational fluid dynamics  CHO  Chinese hamster ovary  CSAF  controlled-shear affinity filtration  CSTR  continuous stirred tank reactor  EBA  expanded bed adsorption  GRM  general rate model of chromatography  HGMS  high-gradient magnetic separation  HIC  hydrophobic interaction chromatography  IEX  ion-exchange  IEXM  ion-exchange membrane  mAB  monoclonal antibody  ODE  ordinary differential equation  PFR  plug flow reactor  [-] [m3 s-1] [rpm]    xxi  PIV  particle image velocimetry  RPLC  reversed phase liquid chromatography  SEC  size exclusion chromatography  SMA  steric mass action  TMP  transmembrane pressure  t-PA  tissue-type plasminogen activator  UF-DF  ultrafiltration-diafiltration  ZRM  zonal rate model  nr-ZRM  non-retention form of the zonal rate model  xxii  Acknowledgements Firstly, I would like to express my gratitude to my supervisor Dr. Charles Haynes for giving me the opportunity to pursue my PhD in a supportive and innovative research environment. His guidance, expertise and encouragement were invaluable not only in the completion of this thesis but also in my development as a scientist. Similarly, I would like to thank my supervisory committee members Drs. Mark Martinez, Fariborz Taghipour and Bruce Bowen for their generosity of time and effort. Their ideas and insights were instrumental in the success of my research. I would also like to give additionally thanks to Dr. Taghipour for offering the use of his laboratory equipment for my particle image velocimetry experiments. Special thanks go to Dr. Eric von Lieres. Without his tireless work and boundless enthusiasm much of my work would have been impossible. In addition I thank him for allowing me to visit his lab in Germany. I would also like to thank the many members, both past and present, of the Michael Smith Laboratories and in particular the Haynes Lab not only for their technical expertise and assistance but also for fostering a friendly working environment. I’d especially like to thank Dr. Louise Creagh for her advice and support. Acknowledgement must also be made to the funding sources that made my work possible, foremost the National Science and Engineering Research Council of Canada. Lastly, I am indebted to my friends (especially those at UBC) for getting me through my PhD with some remaining sanity and to my family for their unwavering support during what has, to them, seemed like my never-ending academic career.  xxiii  Chapter 1  1 Introduction 1.1  Thesis overview  The modern biotechnology era began in the 1970’s with the development of recombinant DNA technology by Boyer and Cohen [1] and hybridoma technology by Kohler and Milstein [2]. Together, these technologies provided a rapid and efficient method for creating microorganisms with specific genetic attributes enabling, among other things, biosynthesis of novel products. The first industrial use of recombinant DNA technology was the production of human insulin, which was approved for medicinal use in 1982 [3]. Since that time more than 200 recombinant biomolecules have been approved for therapeutic use in at least one region of the world [4]. These include blood factors such as Factor IX and VIII for the treatment of haemophilia, hormones such as insulin and human growth hormone, and protein and oligonucleotide based vaccines for the treatment of a wide variety of ailments including hepatitis B, diphtheria and tetanus. Monoclonal antibodies, cytokines such as erythropoietin and interferon, and therapeutic enzymes such as -galactosidase used in the treatment of Fabry disease are also produced by recombinant DNA technology. The global market for these biopharmaceutical products was estimated to be $99 billion US in 2009 [4]. The industrial production of biopharmaceuticals can generally be divided into upstream processing, which includes the culturing of a desired microorganism and the associated production of a target compound, and downstream processing, which consists of the purification of a target compound from either the cell culture medium in the case of a secreted product or from the cell mass itself in the case of an intra-cellular product. Many host cell strains exist for biopharmaceutical production, each offering distinctive cell growth characteristics, as well as product processing and expression capabilities. The choice among bacterial, yeast, plant, insect and mammalian cells is dependent on a number of factors including intellectual property position, the nature and intended use of  1  the target protein, the amount needed, regulatory issues and economic considerations. Of particular note is the requirement in many recombinant protein therapeutics for complex post-translational modifications. One such modification, glycosylation, is present in over 40% of all approved products and affects the protein’s stability, immunogenicity and efficacy [4]. Proper glycosylation and other post-translational modifications are often difficult or impossible to achieve in bacterial, fungal or insect cell lines. Consequently, mammalian cell lines such as Chinese Hamster Ovary (CHO) cells or Baby Hamster Kidney (BHK) cells have become a leading vehicle for biopharmaceutical production, with approximately 70% of all protein pharmaceuticals now produced in these and other mammalian cells [5]. Current technologies for commercial-scale culturing of mammalian cells include roller bottles, microcarriers, and hollow fibre bioreactors cultures [6-8]. However, for largescale manufacturing of biotherapeutics, mammalian cells are now primarily grown in suspension cultures due to the ease of scale-up and more homogeneous and easily controlled properties [6,9]. The simplest growth environment for mammalian cells is provided in batch culture, which involves no feed to or withdrawal from the system during operation. Both cell density and product yields are typically limited in batch cultures by the accumulation of metabolic byproducts, the depletion of nutrients, and/or the inhibition of production by nutrients, metabolites or the product itself. Fed-batch cultures overcome some of these limitations with the addition of nutrient feed over the course of the culture.  This can enhance product yields while retaining operational  simplicity and reliability. Growth and production limitations can also be ameliorated in perfusion cultures, where constant nutrient feed and spent media removal enables a continuous and often high-titre product stream over weeks or months, in part through the removal of labile products from the culture environment.  Though operationally  somewhat more challenging, perfusion bioreactors may be as little as one tenth the size of those required for batch operation [10]. Following their upstream processing, most therapeutic products must be isolated to purity levels above 99%. The downstream processing of protein products from mammalian cell culture typically requires a large number – often more than ten - of sequential unit  2  operations (Figure 1.1). For a secreted protein product, the first stage of purification is cell removal. This is generally followed by a concentration step, usually ultrafiltration (UF), to reduce handling volumes for further purification.  A capture stage is then  employed to separate the target protein from major host and process-related impurities such as proteins, lipids, carbohydrates, proteases, glycosidases, adventitious agents and fermentation byproducts.  Intermediate purification then utilizes high-resolution  techniques such as chromatography to remove most remaining host-cell proteins, endotoxins and nucleic acids. A final polishing stage removes viruses and any closely related products such as degraded protein molecules, while also preparing the product for therapeutic formulation. Production costs of human therapeutic proteins can reach up to $1000 US per gram and, in contrast to early recombinant protein products that were administered in sub-milligram  Step Yield  Overall Yield  Cell culture  100%  100%  Primary Isolation  85%  85%  UF-DF  92%  78%  Capture  90%  70%  Int. Purification  92%  UF-DF Polishing  Figure 1.1  Step Yield  100%  100%  Cell culture  86%  86%  Integrated Isolation & Capture  63%  77%  90%  Int. Purification  92%  58%  71%  92%  UF-DF  95%  55%  67%  95%  Polishing  Comparison of a conventional protein purification stream (left side) to a corresponding bioprocessing scheme where the primary isolation and product capture stages have been integrated (right side).  3  dosages, many current products require dosages of tens or hundreds of milligrams [11,12]. These high costs make some products commercially unfeasible and are the motivation behind efforts to reduce production costs by factors up to 1000-fold. Dramatic improvements in strain engineering, as well as in clone selection, media composition and culturing conditions have pushed typical recombinant product titres in both fed-batch and perfusion cultures as high as 6 g L-1 [5,13,14], and titres of up to 10 g L-1 may become the norm in the near future [11]. Although specific costs of upstream processing are inversely proportional to product titres, these economies of scale are not as pronounced downstream, where costs generally depend more strongly on mass of product than on volume [11]. Consequently, total manufacturing costs for cell-culture derived products are now typically dominated by the relatively high cost of product capture, purification, and formulation, which can account for up to 80% of total production costs [15]. Since downstream processing is the major cost factor, it offers the largest potential for process optimization.  Improvements in downstream processing will have a  significant impact on process efficiency and cost and will offer a competitive advantage in areas where multiple products are marketed for the same indication. The requirement for many downstream processing steps, each of which results in a loss of product, necessarily reduces overall yields. One obvious method to improve product yields is, therefore, to reduce the number of required operations through rational downstream process integration and optimization.  As can be seen in Figure 1.1,  integration of the cell-retention/removal step can potentially increase overall product yield by 10-20% [16]. Integration can also reduce processing times, as well as both operating and capital costs. A promising new technology for the integration of cell separation and product capture is controlled-shear affinity filtration (CSAF).  Preliminary studies on a prototype CSAF  device showed the technology to be effective at capturing recombinant protein directly from mammalian cell culture [16].  Results from this proof of concept study were  obtained without a fundamental understanding of the hydrodynamic, mass transfer and adsorption phenomena within the device, suggesting that significant improvement of the performance of the technology is possible.  Moreover, a sound fundamental  4  understanding of the properties and performance of CSAF is essential for the optimization, scale-up and operation of the technology within a manufacturing pipeline. It is therefore the goal of this thesis to fully develop the CSAF technology by constructing and validating mathematical models that comprehensively describe the hydrodynamics, mass transfer and protein binding kinetics within the device. These models will then represent a theoretical framework that can be used, in the future, in conjunction with experimentally derived data to optimize the design and operation of CSAF for use in an industrial biopharmaceutical production application.  1.2  Current techniques in downstream processing of biopharmaceuticals  Downstream processing of biopharmaceuticals generally consists of several orthogonal unit operations that, together, achieve the rigorous purity standards required of human therapeutics.  Since  most  protein  purification  techniques,  such  as  column  chromatography, cannot tolerate the presence of host cells, downstream processing of cultures of mammalian cells - which secrete product proteins and thus do not require cell disruption – typically begins with a cell separation stage (Figure 1.1). 1.2.1 Cell separation techniques Cell separation can either be designed to recover a clarified supernatant from batch or fed-batch cultures or to continuously retain cells in a perfusion culture. In any case the aim is to deliver a cell-free eluent.  Common cell separation techniques include  sedimentation [17], centrifugation [9,18,19], acoustic filtration [20,21], spin filtration [22-24], flocculation [25] and dynamic microfiltration techniques such as cross-flow filtration [26-29]. Those most commonly used by the biotechnology industry include sedimentation, centrifugation and various modes of microfiltration, as each generally provides good separation efficiency, reasonably long operability, dependable scale-up, and low power requirements.  5  1.2.1.1 Sedimentation Cell retention in gravity settlers operates on the basis of the density difference between cells and the liquid medium. Settlers have been used both for cell removal in batch and fed-batch systems and as retention devices in animal [30], yeast and plant cell perfusion cultures [31,32]. They can be configured as a vertical settling zone [33-35], a single inclined plate [17,36], multiple inclined plates [37-39] or a horizontal settling zone [35]; however, the best results have been obtained in inclined settlers as these devices decrease the sedimentation path and effective upward velocity. Although sedimentation is a simple, robust process that avoids both high shears and the potential for filter clogging, it is limited by the slow settling velocity of mammalian cells caused by their small size (ca. 10 m) and low density (typically only 5% greater than the medium). As a result, high residence times are required, often in excess of 1.5 h [36], and this can result in extended exposure to un-oxygenated environments and harmful proteases that serve to degrade the target biomolecule. Moreover, the separation in gravity settlers is incomplete, producing a clarified supernatant that is not 100% particle free. 1.2.1.2 Centrifugation By applying a centripetal force to increase settling velocities, centrifugation offers significantly shorter residence times compared to sedimentation. Centrifugation systems are closed, compact, scalable, and can be easily cleaned in place. They also can handle a wide variety of feedstocks used in the biotechnology industry. As a result, centrifuges, particularly disk stack centrifuges [40], are widely used for cell harvesting or removal from batch and fed-batch cell cultures.  Centrifugation has also been shown to be  effective in retaining cells in perfusion cultures without decreasing cell viability or growth rate [41,42]. However, mammalian cells can be quite sensitive to shear stress [43-49], making centrifugation unattractive in certain manufacturing operations; for example, some perfusion cultures using a centrifuge as the cell retention device have shown reduced productivity [9,19]. Moreover, the equipment represents a relatively expensive capital investment. Finally, like sedimentation, the harvest is not particle-free,  6  especially when particles less than 1 m are present (i.e. cell debris) [50]. Centrifugation is therefore often operated in conjunction with another clarifying unit operation, such as depth filtration, to produce a 100% cell free eluent [51]. 1.2.1.3 Dynamic microfiltration The most common method for primary recovery of biologics is microfiltration [52]. Microfiltration offers a number of advantages over other cell removal/retention techniques. For one, the membrane barrier in the various modes of microfiltration offers the distinct advantage of providing a particle-free permeate. Microfiltration thus allows for a more consistent feed to the subsequent downstream processing stages, which helps avoid product losses and reduced downstream efficiencies [50]. Furthermore, the high shear environments encountered in centrifugation are avoided in microfiltration and the technology can readily be scaled up from laboratory to industrial-scale volumes. Various configurations and modes of microfiltration are available. The most efficient create fluid motion at the membrane surface to generate shear stress that reduces membrane fouling and the associated need to replace costly membrane cartridges. The simplest method of producing fluid motion at the membrane surface is via forced tangential flow of the cell suspension using either a plate and frame assembly [26,53] or within hollow fibres [28]. However, the generation of this shear stress is linked to pressure gradients within the flow channel that cause inhomogeneous filtrate flux as shown in Figure 1.2 [28]. Several technologies have been proposed to eliminate this problem.  The oldest and most often used is spin filtration [19].  However, that  technology is known to suffer from low cell retention efficiencies [24]. A newer and more elegant strategy is the CSAF technology proposed by Vogel and coworkers [16, 54], which uses a rotating conically shaped rotor to decouple the applied wall shear stress and trans-membrane pressure (TMP), allowing for a uniform permeate flux across the entire membrane surface.  7  Feed  Retentate  Permeate Figure 1.2  Permeate flux characteristics of conventional tangential flow filtration. The transmembrane pressure decreases with position resulting in an inhomogeneous permeate flux.  1.2.2 Protein purification techniques Many factors must be considered when designing a protein purification process. Foremost is the need to maintain the native conformation and biological activity of the target protein while avoiding production of chemical variants through deamidation, oxidation, aggregation or proteolysis. This generally necessitates operating all stages of the process in an aqueous solvent environment and in a manner that avoids extreme temperatures, shears and/or pH levels.  Within these mild and narrow operating  conditions, the purification process must fully remove a wide range of contaminants including host DNA, viral elements, endotoxins and other cellular material. Both the cost and the complexity of a protein purification process are highly dependent on the host cell and cell culture conditions, the size and chemical complexity of the protein product, and the intended use of the product. Although such techniques as twophase aqueous systems [55-59], affinity precipitation [60,61], preparative electrophoresis 62,638], and isoelectric focusing [64] have been used on occasion, purification is most often achieved using a series of classic chromatography steps [65]. As product titres rise in mammalian cell cultures, considerable attention is therefore being given to addressing the well-documented short-comings of classical, bead based chromatography, including  8  through the development and implementation of alternative modes of chromatography such as membrane chromatography [66-71]. 1.2.2.1 Classic packed-bed chromatography Packed-bed chromatography has remained a cornerstone of industrial and therapeutic protein purification for over a half century due to its good capacity, outstanding resolving capability, robustness and reliability. In chromatographic processes the solute of interest is purified from a mixture based on the differential partitioning of the component molecules between the mobile and stationary phases (Figure 1.3). The target molecule can either be retained by the stationary phase while unwanted contaminants are eluted or, conversely, it may be directly eluted from the column as contaminants are retained. Chromatographic stationary phases are distinguished by the functional groups attached to their base matrix. Chromatographic modes that are typically used for the purification of biopharmaceuticals include:   affinity chromatography: relies on specific interactions between the immobilized ligand and its cognate binding partner on the target molecule. Typically offering the highest specificity among available chromatography modes, affinity chromatography also has a high cost. It is used most often for initial product capture [72].    ion-exchange (IEX) chromatography: separates biomolecules based on their unique surface charges at a given pH and ionic strength.  IEX  possesses system robustness, scalability, flexibility, high resolving power and capacity, and is thus widely used for initial capture, intermediate purification and final polishing [72].   hydrophobic  interaction  chromatography  (HIC):  separates  biomolecules based on the interactions between hydrophobic groups on the surface of biomolecules and hydrophobic ligands displayed on the stationary phase. HIC is frequently employed as an orthogonal method to IEX for intermediate purification and final polishing [72].  9  Solute  Increasing Time  Time Figure 1.3  Separation of three solute molecules in elution chromatography. As the mobile phase flows through the column from left to right the three molecules (represented as squares, circles and diamonds) are retained differentially by the stationary phase and thus exit the column at different times.    reversed phase chromatography (RPLC): relies on the partitioning of the solute between the mobile phase and a hydrophobic stationary phase [73]. Although RPLC has a very high resolution, the requirement of organic solvents, which can be denaturing, typically restricts its use to the polishing stages of low molecular weight products [74].    size exclusion chromatography (SEC; gel filtration): separates biomolecules based upon a molecular sieving effect wherein smaller molecules are retained in the porous stationary phase while larger molecules are not and elute quickly. SEC can offer high resolution, but is  10  typically applied only during final polishing because of its lower loading capacity and restricted throughput due to compressible media [74]. An ideal chromatographic separation is characterized by plug flow of the mobile phase, negligible resistances to mass transfer in the mobile and stationary phases, and rapid adsorption and desorption kinetics so that the sorbate is everywhere in equilibrium with free protein in the mobile phase. Consequently, in the case of frontal loading, the protein breakthrough curve (BTC) would be identical in shape to the column input (typically a step change), but shifted in time to account for the loading capacity and time lag due to the finite column volume. An ideal column is therefore infinitely efficient and can achieve the desired separation with no loss of product and maximum loading of the stationary phase. In real chromatographic separations, however, the column efficiency is finite and the BTC is broader in shape compared to the input due to a combination of nonidealities that include: i) axial dispersion effects arising from molecular diffusion and inhomogeneous flow patterns in a packed bed; ii) resistances to solute mass transfer within the mobile phase and to the sorbent surface which may, depending on the configuration of the stationary phase, include diffusion through the fluid film surrounding the stationary phase and within the stagnant fluid in the pore volumes; and, iii) intrinsic kinetics for protein adsorption and desorption that prevent the realization of local equilibrium within the column. This broadening necessarily results in either a loss of product if loading is allowed to proceed past the point of initial breakthrough or a loss of capacity if loading is halted before saturation has been achieved. A central goal of chromatographic design is therefore to sharpen BTCs by reducing the non-ideal effects that lead to band broadening. In preparative-scale bead-based chromatography, adsorption rates are often limited in columns packed with large stationary phase particles by slow solute diffusion into and within the stationary phase, and by pressure-drop limits or binding kinetics in those packed with smaller particles [70]. These limitations represent, along with high costs, the main disadvantages to the technology - although, for high-volume products, scale-up challenges can also be a major disadvantage – and have motivated the development of alternative purification technologies offering higher rates of mass transfer while  11  maintaining high binding capacities and resolving power.  One such technology is  membrane chromatography, which was first introduced in 1988 by Brandt et al. [75] and is gradually gaining industrial acceptance as an attractive alternative to packed bed chromatography for the capture and downstream processing of protein products. 1.2.2.2 Membrane chromatography In membrane chromatography the interconnected pore architecture of the stacked or monolithic column allows solutes to be convected across binding sites (Figure 1.4). Mass transfer is therefore enhanced, limited primarily by the rate of film diffusion or, most often, the intrinsic binding kinetics of the protein-sorbent complex [66]. Because film  Bulk convection  Bulk convection  Film diffusion  Film diffusion  Pore diffusion  Figure 1.4  Solute transport processes in packed bed chromatography (left side) and membrane chromatography (right side). In packed bed chromatography solute molecules must, in order to reach the adsorption surface on the interior of porous beads, diffuse through a thin film layer surrounding the bead and then the stagnant fluid within the pores. Conversely, in membrane chromatography the solute molecules are convected past the binding sites and the only diffusive process is across the thin film layer surrounding the pore wall.  12  diffusion is generally orders of magnitude faster than pore diffusion, mass transport limitations are significantly reduced in membrane chromatography [66]. The technology therefore offers the well-documented potential [67,71,76] to increase throughput and performance provided the column design meets capacity and resolution requirements. Membrane chromatography columns can be configured as tangential-flow, hollow fibre or dead-end filtration systems [66,68,77,78]. In most commercially available units, the membranes are either axially stacked or spirally wound in order to achieve sufficient capacity [66,78,79]. No matter their configuration, membrane columns generally consist of a stack of multiple membranes. This reduces nonideal transport effects arising from local variations in membrane thickness and porosity and sharpens BTCs [70,80].  All  modalities of classical chromatography (i.e. IEX, HIC, affinity and reversed phase) can be employed in membrane chromatography, with ion-exchange membranes currently being the most widely available and used [69]. A multitude of membrane substrates and coupling chemistries have been developed [69,80,81].  Reported applications of  membrane chromatography cover a wide range of target biomolecules, including enzymes, monoclonal antibodies, serum albumin, serum antibodies, DNA, viruses and endotoxins [69,71]. However, because membrane adsorbers are able to handle large volumes of liquid but generally offer a lower static binding capacity, the technology appears best suited to separations where the target protein is either present at lower concentrations in a large volume of feed or is large, so that the minimization of mass transfer by diffusion provides maximum benefit [69]. Advantages to membrane chromatography in addition to enhanced mass transport therefore include increased dynamic capacity for large proteins, as well as ease of scaleup and a simple and inexpensive production process that enables the use of disposable membrane adsorbers – disposability being one of the most promising paradigms to increase efficiencies in biopharmaceutical production [82,83].  Another inherent  advantage of membranes is their sieving capacity: microfiltration membranes are widely used in cell separation (see above) and ultrafiltration membranes are used in buffer exchange and protein concentration.  This sieving capacity, when coupled to an  adsorptive capability, naturally makes membrane chromatography a promising  13  technology for the integration of cell removal and protein capture. The challenge is to implement membrane chromatography in a manner so as to best exploit its distinct advantages. One such implementation, CSAF, is discussed below. 1.2.3 Integrated cell separation and product capture techniques The integration of the cell removal and first product capture stages of downstream processing provides  an excellent  biopharmaceutical production.  opportunity to  improve overall  yields  in  Such process integration can decrease capital and  operating costs, as well as reduce processing times of cell culture bioreactor volumes that are now on the order of 5,000-20,000 L [84]. It is important to keep in mind, however, that integration is only feasible when the performance of the integrated unit is as good or better than the stand-alone units it replaces. Current techniques for the integration of cell removal and product capture include aqueous two-phase extraction, expanded bed adsorption, high-gradient magnetic separation, and controlled-shear affinity filtration. 1.2.3.1 Aqueous two-phase extraction Extraction in aqueous two-phase (ATP) systems has been shown to be an effective means of purifying proteins directly from a cell-containing culture medium at industrial scales [55-59,85-87]. In ATP processes – which can be single or multi-stage - the target protein is preferentially partitioned into one phase while the biomass solids and contaminants (proteins, lipids, carbohydrates, etc.) partition into the other or within the interface. The ATP system can be composed of either two immiscible polymers/detergents or a polymer in a high concentration salt solution; however, only PEG/salt systems have found any meaningful industrial interest due in part to their lower cost, wider applicability, and potential to deliver the clarified extract as a polymer-free solution [56,57]. The use of thermoseparating polymers such as ethylene oxide-propylene oxide is also being investigated to permit polymer removal and recycling and thereby reduce processing costs and carry-over contamination difficulties associated with more traditional phaseforming reagents [87]. ATP extraction has been applied to a range of systems, including the primary recovery and purification of enzymes, therapeutic proteins and monoclonal antibodies [58]. 14  However, ATP extraction is not used widely in industry because of some associated disadvantages. Attaining high product purities typically requires a high concentration of costly phase-forming polymer [57,86], resulting in a high viscosity that hinders further downstream processing and increases costs. Also, the implementation of ATP extraction in continuous mode is not straightforward. Lastly, the process requires either a settling or centrifugation stage to separate the two phases and it typically fails to produce a harvest that is completely cell-free [59]. 1.2.3.2 Expanded bed adsorption Expanded bed adsorption (EBA) is designed to allow for the direct capture of proteins from a cell culture feed without the need for prior clarification [88-92]. Introducing the liquid feed upwards into a bed of unconstrained adsorbent beads causes the bed to expand, creating spaces between the beads that allow for the passage of particulate matter while reducing adsorbent bed clogging and fouling (Figure 1.5). Product capture is achieved through adsorption on the suspended beads, which are washed following loading to remove particulate matter. The suspended bed is then collapsed and the adsorbed product recovered by either isocratic or gradient elution [90-92]. EBA has been shown to be effective in the processing of many cell culture stocks including E. coli, yeast, hybridoma, myeloma and mammalian cells [93], and the technology has been commercialized by GE Healthcare through their Streamline® product line. The columns are equipped with a specially designed distribution system that promotes plug flow and restricts radial flow that would cause turbulence. Most EBA systems utilize an agarose-based sorbent as the stationary phase because of its high chemical and mechanical stability, as well as relatively high binding capacities [93]. An ion exchange modality has been shown to be effective within EBA for the capture of a variety of biomolecules including BSA, plasmid DNA and annexin V [88]. Protein A based affinity systems have also been implemented [90-92,94,95]. While EBA shows promise as a separation technology, there are some significant and well-documented drawbacks. Adhesive cells are problematic as cellular aggregates can clog the column’s flow distributor or become so large that they do not exit the column  15  Upper adapter  Lower adapter  Flow  Figure 1.5  Flow  Comparison of a conventional packed bed chromatography column (left side) to expanded bed adsorption (right side). In conventional packed bed chromatography the bed is constrained and would become plugged with the introduction of particulate matter. In expanded bed adsorption the upward flow of mobile phase through an unconstrained bed of adsorbent beads causes the bed to expand creating spaces between the beads that allow for the passage of particulate matter.  due to an insufficient linear velocity in the upward flow field. Adhesive cells may also bridge sorbent particles to create large stationary phase aggregates, leading to changes in the column’s hydrodynamic properties. Cells may also inhibit target protein binding in certain situations. For instance, anion-exchange chromatography is strongly affected by the presence of negatively charged cells, and a low binding capacity for some affinity EBA resins has likewise been reported [88]. Furthermore, cultures with a high biomass can have physical properties (e.g. viscosity) that are unsuitable for expanded beds [89]. Indeed, to date the technology has generally exhibited a lack of robustness that is characterized by a tendency for dynamic binding capacities to decrease substantially within a limited number of cycles [96]. 16  1.2.3.3 High-gradient magnetic separation High-gradient magnetic separation (HGMS) employs micron-size magnetic adsorbent particles in a scalable batch adsorption step, followed by high-speed recovery of the loaded particles using a magnetic field. Because the particles in HGMS are nonporous, clogging by particulate matter and intraparticle diffusion resistances are eliminated. The technology therefore appears suitable for protein capture direct from cell culture provided reasonable column capacities can be realized. However, while HGMS has been shown effective in laboratory studies of protein capture from culture feedstocks [97-99], it remains an immature unit operation that requires further work, including the development of robust well-characterized magnetic sorbent particles and sterile recovery equipment. 1.2.3.4 Controlled-shear affinity filtration CSAF combines a specially designed rotating disk with an affinity membrane chromatography column to allow for the integrated capture and purification of a secreted recombinant protein product directly from cell culture (Figure 1.6) [16]. The rotor assembly is designed to produce a shear stress at the membrane surface that inhibits membrane fouling and clogging and allows for a radially uniform filtrate flux that serves to minimize broadening of the BTC within the associated affinity membrane chromatography column [68,78]. The shear stress is controlled solely by the rotational speed of the rotor and not, as in the case of more traditional cross flow filtration methods, by pressure gradients (e.g. Figure 1.2).  Therefore, CSAF combines the known  advantages of membrane chromatography, namely reduced mass transport limitations leading to high throughputs in the membrane column [66,75], with the ability of the rotor assembly to decouple the TMP and shear stress at the membrane surface. This permits facile control of the shear stress to which the sensitive mammalian cells are exposed. CSAF therefore offers the potential to maximize cell viability, reducing cell lysis and minimizing contamination from intracellular products such as proteases and DNA. Furthermore, in principal CSAF should have wide applicability.  The technology is  directly amenable to use in batch and fed-batch cultures, and could be used in perfusion systems or any process where a traditional or affinity mode of chromatography can be  17  Figure 1.6  Schematic of CSAF. Feed is introduced to the rotor chamber via A and retentate is removed via B. A conically shaped rotor, E, spins above an affinity membrane stack, D. The shear caused by the rotor inhibits cell deposition allowing for uniform filtrate flux, C.  applied to product capture. However, realizing this potential will require a much better fundamental understanding of the technology than is currently available so as to permit characterization and optimization of the operation for a given application. Although it provides many advantages, the integration of cell removal and product capture as implemented in the original CSAF prototype [16] (Figure 1.6) also poses some potential challenges. Performing cell separation and a chromatographic separation in a fully integrated unit is operationally complex and requires multiple inlets and outlets as well as complicated valving due to the necessity for independently introducing different processing solutions, including cell culture broth, wash buffer, elution buffer and regeneration solution. This implementation is further challenged by the need for tight control of solution conditions (e.g. pH, salt concentration, etc.) in both a bioreactor and the membrane chromatography column.  18  An alternate and arguably more realizable industrial implementation of CSAF is to achieve the desired integration by directly connecting the rotor chamber to a commercially available membrane chromatography capsule as shown in Figure 1.7. The rotor assembly would then be placed above a simple microfiltration membrane for cell retention and removal.  Comprehensive knowledge of hydrodynamics within the rotor  chamber and through the membrane is again required in order to achieve high cell separation performances and a predictable and uniform flow to the membrane chromatography column.  This alternate CSAF configuration loses none of the  advantages of the original prototype technology, but would be much simpler to operate.  B  A  E D C Figure 1.7  Schematic of a connected CSAF configuration in which the rotor, E, spins above a microfiltration membrane. The filtrate from this membrane is then passed directly to an external affinity membrane stack, D. Feed is introduced to the rotor chamber via A, retentate is removed via B and filtrate exits the membrane stack via C.  19  1.3  Thesis objectives  For CSAF to become a feasible alternative to traditional methods of cell separation and product capture, a much better understanding of the hydrodynamic, mass transfer and protein adsorption phenomena occurring in the device is required. Although effective, the prototype CSAF design required operation at high rotational frequencies (ca. 700 rpm) [16]. These high rotational frequencies limit the degree to which the device can be scaled-up and would therefore restrict its use on an industrial scale. It is thus the objective of this thesis to better characterize and refine CSAF technology through the development of mathematical models that comprehensively describe the physical phenomena occurring within the device that influence separation performance. Whichever CSAF configuration is implemented, the technology consists of two interacting parts: a rotor chamber and a membrane chromatography column. Complete characterization of CSAF therefore requires the development of mathematical models that describe: i) the hydrodynamics within the rotor chamber; and, ii) solute mass transfer and protein-ligand binding kinetics within the membrane chromatography column. Here, computational fluid dynamics (CFD) is utilized to model and optimize rotor chamber hydrodynamics. The CFD model discretizes the three-dimensional space and then solves the underlying momentum equations at all positions to give a full picture of the flow fields within the rotor chamber. A new model is then developed to describe solute mass transfer within a stacked membrane chromatography that can either be combined with or directly connected to the rotor chamber. For reasons that will be provided, current models of membrane chromatography often fail to quantitatively reproduce binding and breakthrough data for proteins. Limitations in current modeling approaches are therefore addressed to better predict column breakthrough by specifically accounting for all contributions to system dispersion both within the membrane stack and the extra column hold-up volumes, as well as nonidealities affecting the intrinsic kinetics of the proteinligand binding reaction. In this work the model is validated using experimental data derived from commercially available membrane chromatography modules; however, it is constructed to apply to general configurations of stacked membranes and as such is suitable to describe the membrane stack associated with CSAF.  20  The mathematical models developed in this thesis provide an in-depth description of the fundamental processes occurring during operation of the CSAF technology. Such a description enables the optimization, in silico, of the physical design as well as the definition of the permissible and ideal operational spaces of the device.  These  mathematical models also represent theoretical tools that are essential for technology scale-up, as an ad hoc scale-up methodology consisting of construction and testing of multiple pilot-scale CSAF units is not feasible, at least in a normal academic setting, due to time and cost. Though clearly a long-term goal, the construction and experimental evaluation of a second generation scalable CSAF technology is therefore not the focus of this thesis and would not, in fact, be possible due to prohibitive costs as well as a lack of access to industrial cellstocks. Instead, this thesis lays the theoretical foundation for the future development of CSAF into a useful, industrially applicable alternative to conventional unit operations for the primary recovery and capture of biopharmaceuticals.  1.4  Computational fluid dynamic modeling  Because it can accurately characterize fluid hydrodynamics in complex three-dimensional systems, CFD is finding increasing use in modeling biological processes including bioreactors of various different configurations (e.g. roller bottle [100], tank [101,102], Couette-Taylor [103], membrane [104] and airlift [105,106]) and bioseparation unit operations such as disk-stack centrifuges [107,108], membrane filtration devices [109,110] and chromatography columns [111]. Rotating disk filters have also been studied through CFD modeling [112-115], including recent work by Castilho and Anspach toward modeling a dynamic filter for cell harvesting and recycle that utilizes a rotor assembly above a membrane stack [116]. Their filtration device differs from CSAF in both design and purpose and, as a result, their study did not account for filtrate flux, the presence of an affinity membrane stack, or the need for uniform solvent velocity through the membrane stack. Nevertheless, it demonstrates the potential for using CFD to quantify hydrodynamics in mechanically complex filtration assemblies. In this thesis the commercial software packages FluentTM v. 6.1.22 and GambitTM v. 2.1.6 (Fluent Inc., Lebanon NH) are used for the fluid dynamic calculations and three-  21  dimensional domain discretization respectively. The fundamental mass and momentum conservation equations known as the continuity and Navier-Stokes equations are applied here to an incompressible, Newtonian fluid:    v   0    (1-1)    D  v  P   2 v  g Dt  (1-2)  where, the three-dimensional flow field is given as v ;  and  are the fluid viscosity and density respectively; P is the pressure scalar; t is time; and g represents gravity. These equations are non-linear and cannot be solved analytically except in the simplest of cases. Therefore, Fluent solves these partial differential equations based on the finite-volume method, where the equations are integrated over each control volume yielding discrete equations. Turbulent flow is characterized by small-scale velocity field fluctuations, the exact description of which would require a superfine discretization that would make the model too computationally time-consuming to feasibly solve.  Fluent therefore Reynolds-  averages the instantaneous governing equations by considering solution variables such as velocity and pressure to be the sum of the mean value (designated by capitalization) and the fluctuation (designated by a prime). Thus, for the velocity components:  vi  Vi  vi  (1-3)  Variables of this form are substituted into Equations 1-1 and 1-2 and time-averaged to yield modified forms of the continuity and momentum equations:        V  0 and   v   0  ______    V       V  V   P +  2V    v v  g  t   (1-4)  (1-5)  Equations 1-4 and 1-5 have the same general form as the exact Navier-Stokes equations, with the velocities and other variables now representing time-averaged values. However, 22  they include an additional term denoting the effects of turbulence in the momentum  balance. These additional Reynolds stresses, Rij    vivj  , unlike the time-smoothed viscous stresses, are not related to the velocity gradients in a simple way and are instead complicated functions of position and turbulence intensity. Empirical expressions for the Reynolds stresses are therefore required to solve the governing equations in turbulent flow. Because it is among the simplest and most robust of turbulence models and has been used in similar geometries [115,117], the realizable k- turbulent model is used in this thesis. This model relies on the Boussinesq approach [118,119] wherein the Reynolds stresses are given by:  v v j Rij   t  i   x  j xi   2    k   t vk  3 xk      ij   (1-6)  in which the turbulent viscosity, t, is calculated from:  t   C  k2  (1-7)    In this model the turbulent kinetic energy (k) and its dissipation rate () are determined from the following equations:    k     kvi        t t xi x j  k        vi    t xi x j   t        k     Gk  Gb     YM  S k  x j   (1-8)     2     C S    C  C1 C3 Gb  S  (1-9)  1 2 k k    x j   where,  k 1  v j vi m   , m  S , S  2S ij S ij , S ij    C1  max 0.43,   m  5 2  xi x j        (1-10)  23  Here Gk and Gb represent the generation of turbulent kinetic energy due to the mean velocity gradients and buoyancy, respectively. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, and Sk and S are source terms. C3is a function of the velocity gradients parallel and perpendicular to the gravitational vector and C is a function of velocity gradients and rate-of-rotation tensors. The model constants have been experimentally determined and have been found to work well for a wide variety of flows [118], they are defined as follows: C1 = 1.44, C2 = 1.9, k = 1.0,  = 1.2 The governing equations must also be modified for flow through a porous medium. Because of residence-time requirements flow within the membrane chromatography column occurs at low velocities (~10-4-10-5 m s-1) where laminar flow can be assumed. The membranes are modeled as anisotropic porous media with no inertial losses and a permeability and porosity calculated from data provided by the manufacturer (Pall Inc., New York). Flow through the porous medium is described using the continuity equation and the Brinkman-extended Darcy equation [120]:    u   0  (1-11)        Du   m P   2 u  m u  m Dt K  (1-12)  where m and K are the porosity and permeability of the membranes respectively and u is the three dimensional flow field in the membrane stack. Since CFD results can be highly dependent on the implemented solution strategy (i.e. turbulence model, discretization method, etc.), any CFD model must be evaluated through experiment. The method used in this thesis is particle image velocimetry (PIV). PIV is a powerful non-intrusive technique for obtaining velocities over two-dimensional domains. Pulsing a laser sheet twice in rapid succession illuminates tracer particles within the fluid. This yields two subsequent positions of the particles, which, after postprocessing, provide velocity vectors at any position within the interrogation zone. These  24  results can then be compared to those produced by the CFD model to ensure that the correct form of the governing equations and appropriate boundary conditions are being implemented. The PIV technique has been used extensively as a stand-alone method for fluid flow studies [121-125] as well as in conjunction with CFD models [126-129]. There have also been studies using CFD and PIV to investigate the hydrodynamics in bioreactors [103,106]. The solution scheme of the CFD model is also evaluated at low rotational velocities by comparing model results to a similarity solution derived from von Karman’s work [130] with rotating plates. Such a similarity solution has been used effectively to describe flow in a rotating disk filter similar to the CSAF system [112].  1.5  Modeling of membrane chromatography  Although membrane chromatography can provide enhanced mass transfer compared to traditional packed bed chromatography, it has been used only sparingly in industrial manufacturing settings.  This is because nonidealities caused by axial dispersion,  inhomogeneous radial flow distributions, and relatively high feed and elution manifold dead volumes impact the performance of current membrane chromatography systems [131-133].  Small column heights and low back pressures serve to magnify these  inefficiencies, which must be ameliorated to realize widespread adoption of the technology. As a result, much work over the past two decades has been devoted to developing models for the quantitative description of elution bands and BTCs from membrane chromatography modules [66,70,80,131,133-143]. Several modeling approaches exist for the quantitative description of solute transport and binding within membrane chromatography modules. Most are based on the general rate model (GRM) of nonlinear chromatography which attempts to define the complete set of transport and reaction mechanisms and rates controlling solute breakthrough. In the GRM, solute mass transport within the interstitial volume of a membrane column is described by the one-dimensional continuity equation:  c c  2 c 1   m s  u int  Dx 2  t z  m t z  (1-13)  25  where c(z,t) is the protein concentration in the interstitial volume and depends on axial position (z) and time (t); uint is the interstitial velocity; m is the membrane porosity; Dx is the axial dispersion coefficient; and s(z,t) is the total protein concentration in the stationary phase. This model assumes: i) a constant axial diffusion coefficient, porosity, mobile-phase viscosity and density; ii) non-coupled momentum and mass balance equations (fluid velocity is independent of concentration); iii) mass transport occurs only by axial diffusion and convection and not by other mechanisms such as surface diffusion or radial convection; and, iv) the column is uniformly packed and solvent velocities and solute concentration gradients in the radial direction can be neglected enabling the column to be modeled as a one-dimensional system with a uniform axial velocity profile. Because of the lack of porous particles the GRM can be simplified in membrane chromatography to preclude pore diffusion; however, mass transport resistances due to the stagnant liquid film at the membrane pore surface may be significant and should be taken into account.  The differential mass balance equations to describe such mass  transfer are discussed in Chapter 5. Though relatively comprehensive, the GRM can be operationally cumbersome as a tool for designing and simulating chromatographic separation systems as it requires values for a large number of parameters, some of which are difficult to quantify accurately and unambiguously due in part to differences in the rates of the various transport and binding sub-processes leading to the bound state. However, by accounting for these differences, simplified and typically equally accurate models of chromatography can be derived in which an individual or lumped event is assumed to control the overall rate of protein transport and binding to the stationary phase.  A rigorous mathematical analysis of the  governing equations including the continuity equation and the kinetic binding equations can yield dimensionless groups that provide insight into the relative contributions of the thermodynamic and hydrodynamic effects contributing to BTC broadening.  For  example, Suen and Etzel [70] have previously shown that when the axial Peclet number (Peax = uintL/Dx) is greater than 40, convection dominates solute mass transfer in the interstitial volumes of the membrane stack and the axial dispersion term in the continuity equation can be neglected.  26  Accurate modeling of breakthrough and elution bands for proteins loaded onto membrane chromatography columns also requires an appropriate model describing the intrinsic binding kinetics (and thermodynamics), often termed the binding isotherm model, as well as the accurate determination of all parameters within that model. Selection of an appropriate binding isotherm model is made difficult by the fact that available theories for protein adsorption oversimplify the biophysics of the process, particularly as it occurs in preparative liquid chromatography. For example, the models of binding kinetics and equilibria that are most widely applied to protein adsorption - the individual and multicomponent forms of Langmuir theory - are simple extensions of theories intended to describe adsorption of small non-interacting molecules from the gas phase.  The  adsorption of proteins from aqueous solution onto porous and often chemically heterogeneous solid matrices is considerably more complex. For example, the much greater size and configurational degrees of freedom of proteins, the competition between the protein macro-ion and the solvent for surface binding sites of the same surface energy or differing surface energies, the potential for either orientationally specific or random binding of the sorbate, the ability or inability of the adsorbed protein to diffuse on the sorbent surface, and the non-ideal interactions between components within the solution and bound at the sorbent surface may all contribute to the adsorption rate, energetics and equilibria. A large number of semi-theoretical models for isothermal protein adsorption have therefore been derived to try to account for these non-ideal contributions, including the original and modified forms of the multi-site isotherm [70,131,134], randomsequential-adsorption isotherm [144], spreading isotherm [145,146], and steric massaction (SMA)/hindrance isotherm [147].  The differential equations describing these  models are discussed in Chapter 5. Due to the low length/diameter ratio of stacked membrane chromatography columns, extra-column dispersion effects resulting from mixing, channeling and other processes leading to unequal solute residence times are known to contribute to BTC broadening [132,136,148]. These non-idealities are together referred to as system dispersion and several models have been proposed to account for it, with most employing a plug flow reactor (PFR) in series with a stirred tank reactor (CSTR) and the membrane stack as seen in Figure 1.8 [131,133,142,143,149,150]. Roper and Lightfoot [136] extended this 27  idea with a second PFR/CSTR sequence to also account for dispersion in the eluent collection manifold.  However, it is now known that models that incorporate a linear  PFR/CSTR sequence on either one or both sides of a membrane stack do not in general provide a quantitative description of breakthrough from membrane chromatography modules [133,142,140,149].  The failure of these models arises, in part, from an  oversimplified treatment of transport and binding processes within the membrane stack as evidenced by the recent work of Sarti and coworkers who improved prediction of breakthrough by utilizing a bi-Langmuir isotherm to account for heterogeneities in the binding energy of the sorbent [131]. In addition, although much effort has been made to properly engineer distributors so as to minimize non-idealities, previous studies [132,136,148]  of  flow  within  different  commercially  available  membrane  chromatography units have shown that the average path length through the centerline of the feed-distribution (inlet) manifold is significantly shorter than through outer regions of the manifold, as illustrated in Figure 1.9. Column loading therefore does not occur in the desired plug-flow manner, resulting in solute breakthrough from the centerline of the membrane prior to the saturation of binding sites in the outer radial positions of the membrane stack, even when the axial flow rate is constant and invariant of the radial position.  PFR CSTR Membrane Column Figure 1.8  Diagram of a system dispersion model for stacked membrane chromatography wherein the plug flow reactor (PFR) accounts for a lag time uncoupled from system dispersion and the constant stirred tank reactor (CSTR) describes dispersion effects in the extra column hold up volumes (e.g. tubing, valves, etc.).  28  Q3 Q2 Q1  Q2 Q3  Figure 1.9  Representative path lengths for solute flow within a stacked-membrane chromatography module possessing both axial and radial flow components within the feed distribution and eluent collection manifolds. Each shaded area indicates a radially defined zone in the manifold offering a distinct solute residence time, with the length of the flow line through the zone reflecting the time. Qi represents the flow rate through zone i.  An accurate model-based description of stacked-membrane chromatography requires that all contributions to band broadening be isolated and characterized. In addition to mass transport within the column, the model must therefore also account for all dispersive effects in the extra-column hold up volumes and any non-uniform flow distribution across the membrane stack. By precisely describing all flow and mass transfer related contributions to band broadening, one can then identify a binding model that provides a quantitative description of the protein retention process of interest and determine if that model provides any fundamental insights into the nature of the process, as opposed to simply providing a useful correlation function.  This latter objective could prove  particularly useful to the further development of membrane chromatography of proteins by ion exchange, as important questions remain regarding the relationship between intracolumn events and the elution band profiles typically observed in these preparative separation systems.  In particular, even for well-designed IEX membranes, BTC  asymmetry characterized by a sharp initial breakthrough followed by a slow approach to saturation is often observed [131-133,143,149,150-152]. Moreover, as with traditional packed-bed chromatography, the dynamic binding capacity often depends on linear 29  velocity. However, that dependence is more complex in membrane IEX adsorbers, with a number of studies reporting a non-linear increase in binding capacity with increasing feed velocity within the low to moderate flow regime [149,153,154].  1.6  Thesis outline  This thesis concerns the development of comprehensive mathematical models that describe hydrodynamic, mass transfer and protein adsorption phenomena within both CSAF technology and commercial membrane chromatography systems. These models provide the necessary theoretical tools to optimize the physical design and operating parameters of the CSAF device. They will therefore enable the scale-up of CSAF technology from the prototype level to an industrially relevant scale as well as define the operational space for such a scaled-up unit. As a first step, a CFD model describing the rotor chamber hydrodynamics of a prototypescale CSAF device is developed in Chapter 2, “Optimizing the rotor design for controlled-shear affinity filtration using computational fluid dynamics”. The model is used to refine the design of the rotor geometry, in silico, so as to minimize the operational rotational frequency required to maintain the minimum uniform shear stress level at the membrane surface needed to inhibit cell deposition and achieve uniform filtrate flux through the membrane chromatography column. The CFD model is extended in Chapter 3, “Scale-up of controlled-shear affinity filtration using computational fluid dynamics,” to study the fluid hydrodynamics and expected filter performance of the CSAF device at industrially relevant scales. Additionally, a novel method for the parallelization of CSAF units is presented in this chapter for use in bioprocessing operations larger than 1000 L. A novel zonal rate model (ZRM) for stacked-membrane chromatography is developed in Chapter 4, “Characterizing solute dispersion in stacked membrane chromatography under flow-through conditions using a zonal rate model”. This model improves on existing hold-up volume models to describe dispersion of solute during transport into and out of the column. The model is developed in this chapter under non-binding conditions so as  30  to isolate flow nonidealities associated with the extra-column volumes from any dispersion that is the result of protein binding and desorption. In Chapter 5, “Characterizing stacked membrane chromatography under protein retention conditions using a zonal rate model,” the ZRM is expanded to account for internal mass transfer resistances and intrinsic solute adsorption kinetics and is implemented in a procedural framework that enables the rigorous characterization of all contributions to band broadening. Specifically the ZRM identifies adsorption kinetics as a critical factor effecting column performance and is used to evaluate current kinetic models of solute adsorption, including a novel spreading model.  31  Chapter 2  2 Optimizing the rotor design for controlled-shear affinity filtration using computational fluid dynamics 2.1  Synopsis  The proof of concept study performed by Vogel and co-workers showed that CSAF is an effective method of capturing human tissue-type plasminogen activator (t-PA) and, by extension, other recombinant proteins directly from a CHO culture, producing a 100% cell-free eluate with a product yield of 86% and a purification factor of 16.7 [16]. Inspired by the hydrodynamic properties of the cone-and-plate viscometer, the rotor within the original CSAF device is conical in shape in an attempt to produce a constant and tunable shear stress at the membrane surface that inhibits membrane fouling and clogging by providing both a uniform hydrodynamic force away from [155] and a drag force parallel to [156] the membrane surface. However, at high rotor speeds, centrifugal forces are known to cause an outward radial flow at the rotating surface and an inward radial flow at the plate surface of a cone-andplate rotor geometry [157]. These ‘secondary’ flows are significant when the Reynolds 2  2  [158], where Rrot is the radius of the cone,  is the number, Re, defined as Rrot  rotational frequency of the rotor,  is the angle between the cone and plate, and  is the kinematic viscosity of the fluid, is greater than 1. In the case of the original CSAF device, these radial flow effects are estimated to begin at rotational velocities above 2 rpm and increase with increasing Re. As a result, the original CSAF device, although effective at the 1 L scale, does not provide uniform shear stress at the membrane surface, creating a need for higher than necessary rotor speeds (700 rpm) to achieve the desired performance. Regrettably, high rotor speeds restrict device scale-up, as mammalian cells  32  are compromised by regions of high shear stress and cell lysis will lead to the release of intracellular impurities. The aim of this chapter is to utilize CFD to model the effects of retentate-fluid-chamber and rotor geometries on hydrodynamics in the CSAF device configured in either of the two ways outlined in Section 1.2.3.4. The CFD model of Castilho and Anspach [116] is extended to accurately model hydrodynamics within the fluid chamber and membrane of the CSAF device under conditions where solvent flux through the membrane stack is nonzero. Shear stress profiles at the membrane surface of the CSAF device are computed for a range of rotor geometries constructed in silico to establish an optimal rotor shape. CFD model results are evaluated through comparison with PIV data and results from a simplified numerical solution of the Navier-Stokes equations for conditions of low Re.  2.2  Computational and experimental materials and methods  2.2.1 Computational fluid dynamics The geometry of the CSAF fluid chamber and membrane stack is shown in Figure 1.6, with important dimensions and physical properties of the system provided in Table 2.1. Without affecting the hydrodynamics in the rotor chamber the membrane volume could represent a stacked membrane chromatography column or a microfiltration membrane depending on the desired CSAF configuration.  The partial differential equations  describing mass and momentum conservation within the CSAF fluid chamber are shown in Chapter 1 (Equations 1-1 and 1-2) as are the equations representing flow through the porous medium (Equations 1-11 and 1-12). Equations 1-1, 1-2, 1-11 and 1-12 were solved using a 3-dimensional, double-precision segregated solver (Fluent 6.1.22) in which discretization, pressure-velocity coupling and pressure interpolation were performed using a second-order upwind scheme, the SIMPLE algorithm, and the program PRESTO!, respectively [118,159]. This solution strategy was chosen by considering the geometry of the system, the rotating flow patterns within it, computing time, numerical accuracy and ease of convergence [118,159].  Shear stress  33  Table 2.1  Characteristics of the CSAF device used in the CFD model.  Rotor Chamber Height, Hcham  20 mm  Rotor Chamber Radius, Rcham  43.75 mm  Inlet and Retentate Radius, Rin/ret  3.175 mm  Radial Position of Inlet/Retentate, Rpos  22.9 mm  Rotor Radius, Rrot  35 mm  Rotor Gap Height,   0.2 mm  Membrane Stack Radius, Rmem  30 mm  Membrane Stack Thickness, Ho  450 m  Membrane Permeability, K Membrane Porosity, m  6.6 x 10-15 m2 0.58  profiles at the membrane surface were calculated from the velocity profiles provided by Fluent. The MATLAB code for this post-processing is found in Appendix A. At standard densities of 10 to 20 x 106 cells mL-1, mammalian cells constitute only 2-4% of the total volume of a suspension culture [6]. The working fluid in the simulations was therefore assumed to be homogeneous, Newtonian, isothermal and incompressible. Noslip boundary conditions were applied to all CSAF device walls excluding the membrane surface. For those CFD simulations intended for direct comparison with PIV experiments, where there was no flux through the membrane stack, the fluid chamber was modeled as a closed system with no-slip boundary conditions applied at the membrane surface (retentate side). In all other simulations a 0 Pa outlet pressure was applied at the outlet-flow surface of the membrane stack, and the pressures at the feed inlet and retentate outlet ports were adjusted to achieve the desired filtrate flux. The stackedmembrane column was modeled as a continuous porous medium with an axial viscous resistance two orders of magnitude less than that in the radial and tangential directions. Membrane permeability and porosity were calculated from data provided by the manufacturer (Pall Inc., NY).  All simulations were calculated using a steady-state  assumption in which the membrane permeability remained constant; therefore, membrane  34  fouling was neglected. The realizable k- model was used to describe turbulent flows within the device as this model has been shown to be accurate for similar geometries [115,117]. To ensure the solution method was accurate in fluid regions near moving and non-moving surfaces, an enhanced wall-treatment including pressure gradient effects was implemented.  Iterations continued until all scaled residuals were below 10-4 and  velocities remained constant at representative points below the rotor. The computational grids were built from structured hexahedral elements using Gambit 2.1.6 (Fluent Inc., NH); approximately 700,000 mesh elements were used for each simulation. The grid distribution was densest below the rotor, reaching a maximum just above and within the surface of the membrane (porous region). In the area above the rotor, mesh size was increased with no loss of accuracy.  Mesh independence was  achieved by increasing the number of nodes until no further significant changes were seen in the flow field; increasing the number of mesh elements from 560,000 to 750,000 caused a maximum change in velocity of less than 1%. 2.2.2 Particle image velocimetry PIV measurements were made using a Flow-map 2D system (Dantec Dynamics, NJ) in conjunction with a dual-head Nd:YAG laser (New Wave Research, CA) having a wavelength of 532 nm (green) and energies of up to 52 mJ. The laser sheet was 1 mm thick. A CCD camera (Hamamatsu Photonics, NJ) with a resolution of 1344x1024 pixels and a 12 bit dynamic range was used to capture the images. A Nikon AF Micro-Nikkor lens (60/2.8), fitted with a 514 nm line filter, was used for focusing. The camera and laser were synchronized using the software package Flow Manager (Dantec Dynamics, NJ); the time interval between two pulses was set to 450 s. Each velocity field was calculated from the adaptive cross-correlation of 100 image sets based on interrogation areas of 32x32 pixels with a 50% overlap.  Polyamid seeding particles (Dantec  Dynamics) were used as tracers. These particles had a mean diameter of 20 m and a density of 1.03 g cm-3. The seeding concentration was adjusted so as to achieve between 5 and 10 particles per interrogation area [127].  35  The prototype CSAF device used in the PIV experiments consisted only of the fluid chamber with no membrane stack or associated filtrate flux. It was made of a plexiglass cylinder with an inner diameter of 87.5 mm and a height of 20 mm. The rotor, also made of plexiglass, had a diameter of 70 mm, an angle of 4o, a maximum thickness at the apex of 7.45 mm, and a gap height from the bottom surface of 0.2 mm; it was driven by a variable speed motor. The rotor’s shaft was 5 mm in diameter and was made of stainless steel. In order to decrease distortion of the image caused by refraction of light at curved surfaces, the entire unit was placed in a rectangular plexiglass box filled with glycerine. The laser sheet was orientated vertically to capture tangential and axial velocities along the center-line at a given radial position, or horizontally to capture the associated radial velocity at a given axial position.  2.3  Numerical solution for low Reynolds numbers  In order to evaluate the CFD model, tangential fluid flow profiles within both the fluid chamber and the membrane stack were obtained from an approximate analytical solution of the mass and momentum conservation equations for low Reynolds number flows within a simplified rotor chamber in which the rotor is modeled as having an infinite radius (Figure B.1).  The simplified CSAF-like system therefore consists of a fluid  bounded below by a non-deformable porous membrane and above by an infinite rigid cone of fixed angle of incidence  rotating at a frequency . The distance between the cone and the membrane surface varies with radial position and is defined as h(r). The assumptions made for this model are that the fluid is Newtonian and incompressible, and the flow is laminar; body forces are neglected and the solution is time-independent and axisymmetric. Inertial terms in the equation of motion are retained, as the curvature of the streamlines may be large near the lower boundary. Flow within the membrane stack (porous domain) positioned beneath the rotor was modeled using the Brinkman-extended Darcy equation (Equation 1-12) [120]. The full derivation of this model and its analytical solution for tangential flows can be found in Appendix B. Analytical expressions for axial and radial flow components within both the rotor chamber and the porous medium were not obtained.  36  Results from the model in the limit of no axial flow through the isotropic porous medium (when tangential and radial velocities within the porous medium will be proportionately highest) indicate that, in the case of the CSAF device, tangential flows are negligible at all axial positions of the porous medium. This is confirmed in Figure 2.1, which shows that there are no tangential flows within the porous medium of the simplified CSAF device if the value of Da/m (the Darcy number Da = K/H2, where K is the membrane permeability and H is the total height of the rotor, H = Ho + h(r), divided by the membrane porosity m) is 10-5 or less and fluid flow in the rotor chamber is laminar. Da/m values for the membranes used in the CSAF device are on the order of 1 x 10-9. Furthermore, because the dominant flow in the rotor chamber is in the tangential direction, it can safely be assumed that radial flows in the porous medium are also negligible. CFD simulations performed both in the presence and absence of filtrate flux through the porous medium also show that tangential and radial components of flow are  Figure 2.1  The tangential component of velocity profiles 1 and 2 estimated in the limit of = 0. 1 and 2 represent the dimensionless tangential velocity components in the porous and fluid domains respectively. Four curves are shown: (a) Da/m = 1 x 10-5, (b) Da/m = 1 x 10-3, (c) Da/m = 1 x 10-2 and (d) Da/m = 1 x 105. These estimates were made with B = 0.15,  = 0.25, H = 2 x 10-3 m and m = 0.58.  37  negligible in the porous medium of the CSAF device.  Thus the model derived in  Appendix B could be simplified by replacing the porous medium with a fixed wall allowing a uniform axial filtrate velocity (Um) across its entire surface (Figure 2.2). The numerical solution to this simplified model for all flow components was found by assuming the forms of the velocity components. Following the procedure outlined by von Karman [130] the velocity components may be postulated to take the form:  vr    1 U m    z   r g 2 h(r ) z   h(r )   (2-1)   z   v  r   h(r )   (2-2)   d    z     r 2 z hr   g    z  1 dr r   hr      v z  U m  g  h( r ) 2   h( r )  2       (2-3)  z     h(r) r  Um Figure 2.2  Physical configuration of the 1-domain low Reynolds number numerical solution of the CSAF hydrodynamics. A conical rotor with angle of incidence  and a localized height above the membrane of h(r) is rotated at a frequency . Fluid flows out through the membrane at a constant velocity Um.  38  where g and  are dimensionless functions, to be determined subsequently, that satisfy the continuity equation (Equation 1-1). These equations, after substitution into the momentum equations, elimination of the pressure terms through cross-differentiation of the z- and r-momentum equations and non-dimensionalization, yield the following system of two ordinary differential equations:  d  d  d 2 αh ξ    λη g η  g η λη   2 λη  0  Rrot   dη   dη   dη  (2-4)   d3  4hξ 3 β 2  d   d4 αh ξ   g η 3 g η  λη λη   4 g η  0 3 Rrot αR rot  dη   dη  dη    (2-5)  subject to no-slip and no-penetration boundary conditions at the rotor surface, and no radial or tangential slip at the membrane surface: g 1  g 1  0 g 0   λ0   0  λ1  1 g 0   1  (2-6)  when terms of order h  and smaller have been neglected. Hence this formulation is valid for the case of a nearly constant rotor height over small radial distances. The variables  and  are the non-dimensional coordinates: ξ  r Rrot  (2-7)  η  z hr   (2-8)  and  and  are the Reynolds numbers for flow through the membrane and within the rotor chamber, respectively:  α  U m Rrot    (2-9)  39  β  2 ΩRrot    (2-10)  After specifying a value of the system of two equations was solved using a 4th-order Runge-Kutta shooting method [160]. Trial values of  0 , g 0 and g 0 were assumed and Equations 2-4 and 2-5 were solved as an initial value problem over the domain 0   1. A root finding procedure was then used to update the values of  0 , g 0 and g 0 until the boundary conditions at  = 1 were satisfied to within a  tolerance of 0.1%. MATLAB code for the solution of the one-domain similarity solution can be found in Appendix C.  2.4  Results and discussions  2.4.1 Evaluation of CFD simulations The solution strategy used to model the CSAF device by CFD was evaluated by comparing CFD simulations to results from the numerical model and PIV experiments, respectively. The numerical model is applicable to low rotational velocities corresponding to laminar flow. Past studies have shown for a fixed-angle cone and plate 2  2  ≥ 48 [161]. Note that this Re viscometer that turbulence occurs at Re = Rrot  includes the square of the rotor angle, while  in the similarity solution does not so that it can be generalized to handle non-conical rotor shapes. For a rotor with a fixed incidence angle of 4o, the onset of turbulence would therefore occur at ≥ 9800. For  = 1285 and   = 0.71, where  is the dimensionless rotor radius r/Rrot, Figure 2.3 shows that velocity profiles given by the CFD simulation closely match those of the numerical model. Indeed, good agreement between velocity profiles is observed at all values of  less than 0.75 under laminar flow conditions. Reported shear stress values represent the square root of the sum of the squares of the two shear stress components present at the   membrane surface (  z and  zr ) and were calculated by first fitting the near-wall velocity data to a third-order polynomial using the Levenberg-Marquardt method [162] and then differentiating.  The shear stress profiles calculated from the CFD simulation show good  40  Figure 2.3  CFD computed () dimensionless (A) tangential (), (B) radial (g’) and (C) axial (g) velocity profiles under the rotor at a radial distance from the rotor center-point of  = 0.71, compared to those calculated by the 1domain similarity solution (solid line). Rrot = 35 mm,  = 10 rpm and Um = 3.65 x 10-5 m s-1. 41  agreement with the numerical model (Figure 2.4) at all  up to ca. 0.65, but begin to diverge (as do velocity profiles) as the rotor edge is approached. This divergence is expected since the numerical model applies to a rotor of infinite radius and therefore does not take into account rotor edge effects, while the CFD solution has a finite rotor bounded by a cylindrical wall a small distance away (Figure 1.6). In order to evaluate the CFD model at higher rotational velocities, PIV experiments were performed on a prototype CSAF unit.  As noted before, the prototype CSAF unit  consisted of only the retentate-fluid-chamber and did not include associated fluxes through the membrane cartridge. CFD simulations were therefore carried out for the CSAF device in the absence of filtrate flux, allowing the results to be directly compared with PIV data. In Figure 2.5, representative tangential velocity profiles determined by CFD simulation are compared with PIV data for two different radial positions beneath the rotor. PIV data  Figure 2.4  CFD computed () shear stress profile along the membrane surface compared to that calculated by the 1-domain similarity solution (solid line) when Rrot = 35 mm,  = 10 rpm and Um = 3.65 x 10-5 m s-1.  42  Figure 2.5  Tangential velocity profiles determined by PIV measurement () compared to CFD results (▲). Experiments/simulations were carried out with an angled rotor at a rotational frequency of 700 rpm for radial positions of (A)  = 0.57 and (B)  = 0.86 when  = 0. The error bars represent the standard deviation.  43  are restricted to regions away from the rotor and membrane surface due to the inherent difficulty in capturing velocities near reflective surfaces. Good agreement between the calculated and experimental data sets is observed, which is representative of the level of agreement observed at other radial positions and rotational velocities. Thus, the CFD results are seen to accurately represent both the experimental and modeling data, indicating that the solution strategy was effective and that the realizable k- turbulence model was appropriate over the range of permissible CSAF operating conditions. The CFD calculated tangential velocity profiles indicate that there are two boundary layers in the gap region: one near the rotor surface and the second adjacent to the membrane (stator) surface. In the intermediate region between the two surfaces, the fluid rotates at a nearly constant velocity Xr, where 0<X<1. This type of velocity profile has been reported previously for flow geometries similar to that beneath the CSAF rotor [163, 164]. In addition, although much weaker than tangential flows, a radial outflow of fluid is observed along the rotor surface and an inflow at the stator surface. 2.4.2 CFD modeling of original CSAF device Once evaluated, the CFD model was used to investigate the shear stress profile at the membrane surface of the original CSAF device, which utilizes a rotor having a fixed angle of 4° and a gap height of 0.2 mm between the rotor tip and the membrane surface. When the system is operated at a filtrate flux of 125 L m-2 hr-1 and its optimal rotor speed of 700 rpm [16], CFD simulations show that the shear stress at the membrane surface has a complex dependence on radial position, increasing to a local maximum near  = 0.1 and then declining before increasing abruptly at  above ca. 0.55 (Figure 2.6). The presence of filtrate flux lessens the formation of the boundary layers (Figure 2.7) compared to what was observed in the absence of filtrate flux (Figure 2.5).  As a result, smooth  velocity profiles are observed under a larger portion of the rotor, with the onset of boundary layer formation at the membrane and rotor surfaces occurring at  ca. 0.55. This emergence of distinct boundary layers, and the associated steepness in the tangential velocity profiles at  ≥ 0.55, coincides with the abrupt increase in shear stress (Figure  44  Figure 2.6  Shear stress profiles along the membrane surface for the original CSAF rotor at  = 700 rpm () and for the CFD optimized rotor at  = 250 rpm (▲) when Rrot = 35 mm and filtrate flux is 125 L m-2 hr-1.  Figure 2.7  Tangential velocity profiles under the original CSAF rotor at various radial positions: (a)  = 0.29, (b)  = 0.57, (c)  = 0.71 and (d)  = 0.86 when  = 700 rpm, Rrot = 35 mm and the filtrate flux is 125 L m-2 hr-1.  45  2.6), indicating that boundary layer formation must be prevented or at least minimized to achieve constant shear stress across the entire membrane surface. In the original CSAF device, Rmem = 0.85Rrot; that is, the radius of the membrane Rmem is less than that of the rotor Rrot to minimize distortions in shear stress profiles at the membrane surface due to rotor edge effects.  Although constant shear stress at the  membrane surface was not achieved, the original CSAF device was shown to be effective at both cell removal and product capture. When combined with the CFD results, this indicates that a threshold shear stress value exists for a given filtrate flux above which cells are effectively swept off of the membrane. At a rotor speed of 700 rpm and a filtrate flux of 125 L m-2 hr-1, the shear stress falls to a local minimum of ca. 0.17 Pa near   = 0.5, suggesting that the satisfactory performance of the original CSAF technology is due to meeting or exceeding this threshold shear stress at all radial positions while at the same time creating a constant pressure at the membrane surface to ensure a homogeneous filtrate flux. The CFD results show that the latter condition is met, as the variations in pressure at the membrane surface are more than two orders of magnitude less than the transmembrane pressure. As the high rotor velocities required to achieve a threshold shear stress of 0.17 Pa adversely affect scalability of the original CSAF system, the use of in silico CFD simulations to reshape the rotor was explored such that constant shear stress at or above 0.17 Pa was observed across the membrane surface at significantly reduced rotor speeds. 2.4.3 Optimization of CSAF rotor geometry The geometry of the CSAF rotor was optimized by adjusting the shape of the original rotor based on the calculated shear stress profile at the membrane surface and a secanttype root-finding method based on the approximate simplifying assumption that the local shear stress at a given radial position r decreases linearly with h(r), the distance between the rotor and membrane at r. The rotor height was thereby decreased by a calculated weighting factor at radial positions characterized by low shear stress at the membrane surface, and vice versa. Throughout the optimization process, the solvent flux through the membrane was set at 125 L m-2 hr-1, consistent with the volumetric flow rate through  46  the membrane during steady-state operation of the original device. The hydrodynamics of the newly designed rotor were then determined by CFD simulation, and the process repeated (7 iterations in total) until a design was achieved that produced a constant shear stress profile everywhere except near r = 0; r is zero at the center-point of the rotor irrespective of . As a result, the shear stress tends to zero as r approaches zero. The shape of the CFD optimized rotor, which is compared to the original rotor design in Figure 2.8, is defined by the equation: h(r) = 1.84·105m-4 · r5 - 2.02·104m-3 · r4 + 8.28·102m-2 · r3 … – 1.57·101m-1 · r2 + 1.53·10-1 · r + 2.00·10-4m  (2-11)  Due to its variable angle of incidence with the membrane, the new rotor creates a significantly more desirable shear stress profile at the membrane surface (Figure 2.6). Furthermore, because of the reduced rotational velocity and the refinement in rotor shape, the onset of boundary layer formation in the tangential velocity profiles (Figure 2.9) occurs at  ca. 0.86; that is, at a radial position beyond the membrane. As a result, no abrupt changes in shear stress are observed. Instead, uniform shear stress is achieved at all radial positions except directly beneath the rotor center-point and at the outer edge of the membrane stack where shear stress increases.  Figure 2.8  More importantly, the threshold shear  Schematic of (a) the CFD optimized rotor shape compared to (b) the original CSAF rotor having a fixed angle of 4o.  47  Figure 2.9  Tangential velocity profiles under the CFD optimized rotor at various radial positions: (a)  = 0.29, (b)  = 0.57, (c)  = 0.71 and (d)  = 0.86 when  = 250 rpm, Rrot = 35 mm and the filtrate flux is 125 L m-2 hr-1.  stress value of 0.17 Pa is achieved across virtually the entire membrane surface at a rotor speed of ca. 250 rpm, 64% less than the rotor speed required in the original CSAF device. Consequently the new rotor should greatly improve the ability to scale-up the technology. In addition, with this change in rotor geometry, the pressure remains essentially constant across the entire membrane surface, varying by less than 1% of the TMP and thereby providing uniform filtrate flow through either the affinity membrane stack or the microfiltration membrane depending on CSAF configuration. It is important to note that the objective function for rotor optimization sought to achieve a uniform surface shear stress of 0.17 Pa at the lowest possible rotor velocity and at a filtrate flux of 125 L m-2 hr-1. As different host cells or culture conditions could require the CSAF to be operated at different rotor speeds or filtrate fluxes, the performance sensitivity of the optimized rotor to these two process variables was investigated. Increases in filtrate flux up to 50 times that used in the original CSAF studies [16] had no effect on the shear stress profile at the membrane surface due to the fact that Um remains much smaller than the tangential velocities in the rotor chamber. The effect of the 48  rotational velocity on the performance of the optimized rotor is shown in Figure 2.10. Rotor velocities between 200 and 500 rpm yield essentially uniform shear stress at the membrane surface, with the value of the surface shear stress increasing with increasing rotor speed. At 700 rpm, the shear stress at the membrane surface is less uniform, showing a small local maximum followed by a shallow local minimum, due to the increased effect of turbulent secondary flows. Thus, the new rotor produces the desired shear stress profile for a wide range of rotational velocities. However, if the rotational velocity required to achieve sufficient surface shear stress to clear retained cells from the membrane lies below 100 rpm or above 500 rpm, different rotor geometries may be required. The CFD simulations reported here provide an efficient and effective means for designing such rotors.  Figure 2.10  Shear stress profiles along the membrane surface for the CFD optimized rotor at different rotor speeds: 200 rpm (), 300 rpm (▲), 400 rpm (), 500 rpm (+) and 700 rpm (). Rrot = 35 mm and filtrate flux is 125 L m-2 hr-1.  49  2.5  Concluding remarks  Close agreement of velocity profiles determined by CFD simulations with those determined by PIV experiments and an appropriate numerical model indicate the capability of CFD for modeling the hydrodynamics within the rotor chamber of the CSAF technology.  CFD simulations of the original CSAF design show that  unnecessarily high rotor velocities are required to achieve sufficient shear stress to clear cells from the membrane surface and ensure a uniform filtrate flux and pressure over the entire membrane surface.  At the bench-top-scale, such high rotor speeds are not  problematic; however, they become disadvantageous during scale-up as they will limit rotor size: the larger the rotor the higher the rotor tip velocity, resulting in regions of high shear stress that promote cell lysis. CFD simulations were therefore used to redesign the rotor to attain comparable and uniform surface shear stresses at lower rotational velocities. The CFD optimized rotor provides uniform surface shear stress, reaching the required threshold value of 0.17 Pa at a rotational velocity that is 64% lower than that required in the original device. The reduced rotational velocity will permit the maximum allowable rotor diameter to be increased by a factor of 2.8 over that possible in the original design. This will result in almost an eight-fold increase in the effective filtration area and, thus, the volumetric throughput of the device. Moreover, the CFD model framework developed herein will enable further in silico design optimization and scaleup and therefore represents an essential theoretical tool that will allow CSAF to become a useful alternative to stand alone cell removal and product capture stages in the industrial purification of recombinant protein therapeutics.  50  Chapter 3  3 Scale-up of controlled-shear affinity filtration using computational fluid dynamics 3.1  Synopsis  In Chapter 2 it was demonstrated through CFD modeling that the scale-up potential of the original CSAF design was limited by unfavourable hydrodynamics that required operation of the device at high rotational frequencies. The CFD model was used to improve the laboratory-scale CSAF design so that minimal cell deposition and desired filtrate fluxes were obtained at a rotational frequency of 250 rpm. This design allowed for a degree of scale-up almost eight times that which could have been obtained in the original device. While clearly an improvement, this still did not represent a scale that would permit the application of CSAF in a manufacturing setting. The aim of this chapter was therefore to refine and expand the original CFD simulations to permit full optimization of the CSAF device to achieve a design that will fully process a 1,000 L bioreactor in a timeframe (ca. 2 hours) competitive with current cell harvesting technologies. CFD was used to construct, in silico, CSAF devices with rotor sizes up to 140 cm in radius. The fluid dynamics within the rotor chamber of larger-scale CSAF units are shown to be complex and include turbulent boundary layers indicating that CFD likely provides the only reliable route to CSAF scale-up. The CFD models were then used to determine the changes in rotor and chamber geometry required during scale-up to achieve the minimum rotational frequency necessary to effectively inhibit cell deposition at the membrane surface. These results then permitted the determination of the maximum rotor diameter possible before reaching rotor edge shear stresses known to reduce mammalian cell viability. This limit was therefore used to set the maximum size and throughput of a single CSAF device. Finally, a configuration for the parallelization of  51  CSAF units is proposed in order to accommodate scale-up beyond the throughput limit of a single device.  3.2  CFD model and computational methods  A schematic for the scaled-up CSAF is seen in Figure 3-1 and important dimensions are listed in Table 3-1. As in Chapter 2, the membrane volume can represent either a stacked membrane chromatography column or a microfiltration membrane depending on the CSAF configuration (see Section 1.2.3.4).  The rotor gap height and rotational and  filtrate velocities are of the same order in the scaled-up CSAF as in the small-scale device. As a result, the CFD model and solution strategy that was evaluated and verified in Chapter 2 can be used as a basis point for the simulations reported here for device optimization and scale-up. This is important, as the original CFD model was chosen based on careful consideration of the geometry of the system, the rotating flow patterns within it, the computing time, the numerical accuracy and the ease of convergence [118,164]. A full description of the CFD model, grid meshing and solution strategy can be found in Section 2.2.1.  Figure 3.1  Schematic of a scale-up appropriate CSAF device. A and B refer to the feed and retentate lines respectively and have been shifted from the top of the unit to the sides. C refers to the filtrate from the membrane stack, D. The rotor, E, is supported by a shaft that now extends through the entire CSAF device.  52  Table 3.1  Characteristics of the scaled CSAF devices used in the CFD model.  Rotor Chamber Height, Hcham  2.0 cm  1st level of Scale-up 20 cm  Rotor Chamber Radius, Rcham  4.375 cm  43.75 cm  180 cm  Shaft Radius, Rshaft  0.25 cm  2.5 cm  27 cm  Rotor Radius, Rrot  3.5 cm  35 cm  140 cm  Max Rotor Height, hmax  0.11 cm  0.40 cm  0.60 cm  Membrane Stack Radius, Rmem  0.30 cm  30 cm  140 cm  Membrane Stack Thickness, Ho  0.045 cm  0.45 cm  3 cm  6.6 x 10-15 m2  6.6 x 10-15 m2  6.6 x 10-15 m2  0.58  0.58  0.58  Small-scale  Membrane Permeability, K Membrane Porosity, m  2nd level of Scale-up 80 cm  In the scaled-up CSAF device fully turbulent regions exist where the turbulent Reynolds number, Rey, was above 200.  In such regions shear stress calculations used the turbulent  viscosity [118]. In regions where the near wall mesh was sufficiently fine to resolve the viscous sublayer (i.e. where the dimensionless wall distance, y+ ~ 1) and where Rey was below 200, shear stresses were calculated using the laminar viscosity.  3.3  Results and discussions  As shown in Chapter 2 CFD modeling of the original CSAF device resulted in refinements to the rotor and chamber design to create a second-generation CSAF unit in which the shear stress at most positions of the membrane surface (Figure 3.2) was uniform and above 0.17 Pa for rotational frequencies () of 250 rpm or greater. This shear stress threshold of 0.17 Pa is important, as it was previously determined to be sufficient to maintain a cell-free membrane surface at linear velocities through the affinity membrane column up to 0.12 cm min-1 [16], which encompasses the range of linear velocities over which the CSAF unit typically operates. However, the shear stress in these second-generation CSAF units is not uniform over the entire membrane surface (Figure 3.2). At the outer edge of the membrane stack, edge effects cause the shear stress  53  Figure 3.2  Shear stress profile along the membrane surface for the optimized, smallscale CSAF device.  = 250 rpm, Rrot = 3.5 cm.  to increase. This is not problematic as the shear stress is greater than that required to maintain a cell-free membrane surface. What is problematic is the rapid decrease in the shear stress at the membrane surface beneath the rotor centerline, which is due to the drop in fluid velocity beneath this slower moving (since vrot(r) = r) region of the rotor. This results in an area near the centre of the membrane stack where cell deposition occurs and a significant decrease in the linear velocity of filtrate through the membrane is observed. This in turn leads to a diminished dynamic capacity and significant broadening of product breakthrough. To overcome this problem, and as a first step towards scale-up, the CFD model was used to simulate the result of extending the rotor shaft in the second-generation CSAF unit through the membrane stack, as seen in Figure 3.1. Although the implementation of a ring-shaped membrane stack surrounding a non-permeable inner core would reduce effective filtration area, the extension of the rotor shaft provides certain benefits.  The  primary advantage of this proposed third-generation design is that it removes the centerline portion of the membrane where cell deposition occurred. Furthermore, by  54  coupling this change in rotor shaft design with a relocation of the inlet and retentate lines from the top of the unit to the sides (Figure 3.1), the CSAF units can be parallelized in a novel manner (see below) to facilitate scale-up. CFD simulations were used to refine the rotor shape for this new shaft geometry and for the accompanying shift in the position of the inlet (feed) and outlet (retentate) ports. Focusing first on a laboratory-scale version of the device (Rrot = 3.5 cm), the goal of the simulations was to identify a rotor contour and positioning that minimizes the rotation speed required to achieve a uniform shear stress greater than 0.17 Pa over the entire membrane surface. As shown in Figure 3.3, CFD optimization of the updated design results in a shear stress profile at the membrane surface that no longer drops to zero at the centerline. Rather, as with the outer edge of the rotor, the shear stress increases at membrane positions adjacent to the rotor shaft, suggesting that all previous areas of cell deposition have been removed in the new design. Moreover, by properly contouring the lower surface of the rotor, a constant shear stress across the entire surface of the affinity membrane column is achieved.  Figure 3.3  Shear stress profile along the membrane surface for a small-scale CSAF device that has side feed and retentate inlets and a shaft extending through the entire device.  = 250 rpm, Rrot = 3.5 cm and Rshaft = 0.25 cm.  55  These encouraging design changes to the laboratory-scale CSAF unit were then used as a starting point for applying the CFD model to scale-up the CSAF device in silico to a 10fold increase in Rrot to 35 cm. All other CSAF dimensions were scaled similarly (Table 3.1). The increase in rotor size was not coupled to a proportional decrease in rotational frequency, thus resulting in a 100-fold increase in the rotational Reynolds number, . Consequently, the shear stress profiles at large-scale are not expected to match those from the small-scale device. In fact, exhaustive CFD simulations across a wide range of geometries of the Rrot = 35 cm CSAF unit revealed that no combination of rotor shape and position results in a uniform shear stress across the entire membrane surface. As a result, the previously described methods used in Chapter 2 to generate a constant surface shear in laboratory-scale CSAF systems do not apply towards design of large-scale CSAF units. Instead, while a relatively small area of constant surface shear stress can be created across the membrane region beneath and close to the rotor centerline (Figure 3.4 inset), a second hydrodynamic domain where the surface shear stress increases with radial position r is always observed. A typical example of this phenomenon is shown in Figure 3.4 for a  of 250 rpm, where the constant shear stress "inner" domain encompasses the portion of the membrane surface extending from the rotor shaft to  = ca. 0.15. The second "outer" domain stretches from  = ca. 0.2 to the outer edge of the membrane stack, within which a nearly linear increase in shear stress with r is predicted. The CFD simulations reveal that the appearance of this second domain in larger CSAF devices is due to the formation of turbulent boundary layers and that the transition from the inner to the outer domain is always characterized by an abrupt increase in shear stress. The inner domain is characterized by a laminar sublayer in which tangential velocities (the major component of shear stress) increase gradually and monotonically with distance from the membrane surface (Figure 3.5). This leads to relatively low shear stresses at the membrane surface that can be, as shown previously, made independent of radial position through proper contouring of the rotor. In contrast, at radial positions greater than  = 0.2, tangential velocities increase sharply with axial position away from the membrane and rotor surfaces due to the formation of a turbulent boundary layer at each surface  56  Figure 3.4  Shear stress profile along the membrane surface for a CSAF device that has been scaled-up by a factor of 10 fold.  = 250 rpm, Rrot = 35 cm and Rshaft = 2.5 cm. The inset shows that the minimum shear threshold of 0.17 Pa (dotted line) is being maintained over an inner domain that stretches from the shaft to a radial position of  = 0.15.  Figure 3.5  Dimensionless tangential velocity profiles under a scaled-up CSAF rotor at various radial positions: (a)  = 0.11, (b)  = 0.43 and (c)  = 0.80 when  = 250 rpm and Rrot = 35 cm.  57  (Figure 3.5). The velocity gradient in these turbulent boundary layers depends on radial position, resulting in a radial dependence to the shear stress at the membrane surface. This dependence, as well as the possible presence of two distinct hydrodynamic domains, complicate scale-up, as they necessitate determination of radial, tangential and axial velocity vectors at all positions between the rotor and membrane surface, as well as their dependence on the geometry of the rotor-containing feed chamber (e.g. rotor diameter and contour, shaft diameter, , etc.). Experimental determination of this information is at best difficult. But it can be obtained through CFD modeling, making the in silico testing of putative CSAF designs by CFD simulation a logical first step toward technology scale-up. It was therefore sought to use CFD simulations to first determine the dependence of the inner and outer domain hydrodynamics on CSAF design, scale and operating conditions, and to then define the changes in CSAF design required for system scale-up. As noted above, in larger CSAF devices no combination of rotor shape, rotor positioning and  eliminates the outer hydrodynamic domain to provide a constant shear stress across the membrane surface.  Scale-up of CSAF therefore requires knowledge of  velocity profiles and shear stresses in the outer hydrodynamic domain. It is however important to recognize that a rapid increase in surface shear stress is always observed in the transition region between the two hydrodynamic domains (Figure 3.4), which causes the lowest surface shear stress in the outer domain to be significantly greater than the uniform shear stress that can be created within the inner hydrodynamic domain.  The  exploration of CSAF performance by CFD simulations therefore provided an intriguing new avenue for scale-up based on searching for designs that eliminate the inner hydrodynamic domain. The threshold shear stress requirement of 0.17 Pa can then be achieved at much lower  values, allowing one to increase the rotor radius, and thus the total cross-sectional area of the affinity column, until the point where the shear stress created at the outer edge of the rotor reaches a set point below that which is known to diminish cell viability. The level of acceptable shear stress for mammalian cells has previously been shown to depend on cell line as well as the type and duration of the applied shear stress [43-47]. However, taken together, the data indicate that mammalian  58  cell viability is unlikely to be compromised by shear stress values less than 70 to 80 Pa, particularly in cases where exposure to such stress is rare and short-lived, as is the case in the CSAF system [48,49]. To establish a conservative design, the maximum allowable shear stress in a CSAF device is therefore set as 0.75 x 70 Pa (= 52 Pa). To illustrate the potential value of the new CSAF design concept, the shear stress profile shown in Figure 3.4 for operation of the unoptimized third-generation CSAF device (Rrot = 35 cm) is first considered. A high rotational frequency of 250 rpm is needed to generate a uniform threshold shear stress of 0.17 Pa within the inner hydrodynamic region extending from the shaft radius ( = 0.07) to the inner/outer domain transition ( = 0.15).  This high  value severely limits scale-up of this CSAF device since the  maximum possible Rrot for this rotor geometry, set by the maximum allowable shear stress, is a mere 11 cm (equivalent to  = 0.31 in Figure 3.4). Now imagine that the inner hydrodynamic zone can be eliminated from the device by increasing the shaft diameter of the rotor to  ~ 0.2. At  = 250 rpm, the minimum shear stress would then be ca. 4 Pa, well above the 0.17 Pa threshold value required to prevent cell deposition on the membrane surface. A significant reduction in  might therefore be possible with this new CSAF design to allow for a significant increase in Rrot and, consequently, an increase in the total cross-sectional area of the affinity column. The goal then was to use CFD to search for combinations of rotor shaft diameter, rotor shape, and rotational speed that satisfy this design requirement and thereby maximize the cross-sectional area and throughput of the affinity membrane column. Simulation results for one such combination where the rotor shape is described by: h(r) = 1.94·10-3m-2 · r3 – 5.71·10-3m-1 · r2 + 6.56·10-3 · r + 2.70·10-3m  (3-1)  are shown in Figure 3.6, where an increase in the shaft radius to 27 cm completely eliminates the inner hydrodynamic domain and thereby allows  to be reduced to 20 rpm while still achieving a minimum surface shear stress that exceeds the 0.17 Pa threshold value for successful CSAF operation. This more than one order of magnitude decrease in  59  Figure 3.6  Shear stress profile along the membrane surface for a CSAF device that has been scaled-up by a factor of 40 fold.  = 20 rpm, Rrot = 140 cm and Rshaft = 27 cm.   permits the rotor radius to be increased to 140 cm without generating shear stresses at the rotor edge that could diminish cell viability. If the outer radius of the membrane stack is aligned with the rotor edge (i.e.  = 1.0), this novel third-generation design would therefore provide an effective filtration area of 5.93 m2 which would represent a more than 1,000-fold improvement in the maximum volumetric throughput of a single CSAF unit when compared against the 0.005 m2 maximum cross-sectional area that can be achieved through scale-up of the original prototype design operating at 700 rpm. While the minimum shear stress threshold is comfortably exceeded across the entire membrane surface (Figure 3.6), the true feasibility of this third-generation design depends on a more comprehensive understanding of the non-uniform shear stress profile created at the membrane surface. Of particular concern is the need to generate a uniform pressure field across the entire membrane surface so that linear (axial) fluid velocities are constant at every radial position within the affinity membrane column. As the computed shear stress at each position on the membrane surface (Figure 3.6) is a component of the absolute pressure at that position, the dependence of the shear stress on r could result in a  60  non-uniform absolute pressure acting on fluids entering the membrane stack. However, the shear stress created by the action of the rotor is quite small compared to the applied pressure field used to pump feed into the rotor chamber. As a result, absolute pressures at the membrane surface vary by less than half a percent (Figure 3.7). Axial fluid velocities through the affinity membrane column are therefore uniform (Figure 3.8), creating favorable plug flow conditions in the column that should lead to efficient product capture and sharp breakthrough. Although membrane fouling and associated loss of permeability are always concerns in manufacturing, past results at the prototype-scale indicate that those issues should not be problematic for the CSAF system. Vogel et al. demonstrated that TMP increases only slightly across twelve successive cycles [16].  Furthermore, this increase is the result not  of membrane fouling but of the increase in cell density in the rotor chamber due to the concentration of cell containing culture medium as it is processed by the CSAF device. Since the shear stresses at the membrane surface in the scaled-up device are higher than  Figure 3.7  Normalized absolute pressure profile along the membrane surface for a CSAF device that has been scaled-up by a factor of 40 fold.  = 20 rpm, Rrot = 140 cm and Rshaft = 27 cm.  = P(r)/Pmax where P(r) is the absolute pressure and Pmax is the maximum absolute pressure seen at the membrane surface.  61  Figure 3.8  Normalized axial velocity profile in the affinity membrane stack for a CSAF device that has been scaled-up by a factor of 40 fold.  = 20 rpm, Rrot = 140 cm and Rshaft = 27 cm. ax = uax(r)/uax,max where uax(r) is the axial velocity in the membrane stack and uax max is the maximum axial velocity in the affinity membrane stack.  those in the prototype it is expected that fouling and any associated decrease in membrane permeability will not occur. Experiments conducted with the original CSAF system show that maximum filtrate flux increases approximately linearly with increasing shear stress, with a maximum flux of 70 L m-2 h-1 observed at 0.17 Pa [16]. Taking then a value of 100 L m-2 h-1 as a conservative basis for design at 0.34 Pa (i.e. equivalent to that generated in the large-scale third generation device), the membrane stack (cross-sectional area of 5.93 m2) beneath the 140 cm rotor offers a total cell-free supernatant loading capacity of 590 L h-1. For the combined CSAF configuration (Figure 1.6) a membrane stack between 0.09 and 0.9 cm thick, which is easily realized in CSAF units, would then be capable of fully processing in ca. 2 hours up to 1,000 L of culture with a product titre of 5 g L-1 based on a typical dynamic capacity within adsorptive membrane columns of between 5 and 50 g L-1 [77,78,165]. A 2-hour processing time is comparatively short when compared to typical  62  doubling times for mammalian cell growth. A single CSAF device could therefore be used in a pseudo-continuous manner for harvest in perfusion cultures as well as in smaller volume batch and fed-batch modes of operation. The handling of larger throughputs would, however, require parallelizing CSAF units to avoid the cell-damaging shear stresses that a larger rotor would create. Here again one can take advantage of the CFD simulations and the associated discovery of the novel extended rotor shaft geometry that will serve as the basis for the third-generation technology. By repositioning the feed inlet and the retentate outlet to the side of the rotor chamber, scale-up through parallelization can be achieved by stacking CSAF units one atop another as shown in Figure 3.9. The filtrate collectors of the individual CSAF units – which must be designed to minimize volume so that the eluent is not diluted and the concentration factor remains high – can then either allow for independent processing or pooling all the units’ filtrates together in a combined outflow. There are however, certain benefits to independent processing, namely that a failure of one membrane stack will not affect an entire product batch. The advantages of the stacked configuration are that it has  Figure 3.9  Proposed configuration for the parallelization of CSAF units. Individual CSAF units are stacked one atop the other with one centralized shaft and individual feed (A), retentate (B) and filtrate (C) lines.  63  a low footprint, it requires only one central shaft and therefore one centralized drive system, and finally, it offers, for the combined CSAF configuration, low cycle times because the large total bed volume is broken up between a series of individual lowervolume CSAF units.  3.4  Concluding remarks  Though CSAF technology has already shown promise as a novel integrated technology for cell retention and product capture and purification from mammalian cell culture, its application in industry requires an effective method for technology scale-up. In Chapter 2 comparison with experimental flow and particle image velocimetry data was used to prove the ability to construct in silico CFD models of CSAF that accurately predict hydrodynamics within the device.  Here, these CFD models have been used to show that  the hydrodynamics within the CSAF rotor chamber (Figure 3.1) are extremely complex and exhibit a strong nonlinear dependence on the size of the rotor and flow chamber. This important result, when combined with the ability to use CFD simulations to rapidly screen a wide range of putative CSAF rotor-chamber designs, strongly suggests that CFD modeling represents the most effective first step in scale-up of CSAF technology. Toward this goal, CFD modeling was used in Chapter 2 to discover a design flaw in the original CSAF device that created an undesirable need to operate at a  of 700 rpm in order to maintain a required minimum shear stress threshold of 0.17 Pa at the membrane surface. Based on this requirement and taking into consideration the shear sensitivity of mammalian cells, this initial prototype could therefore have been scaled-up to an effective filtration area of only 0.005 m2. In this Chapter, CFD modeling has identified the formation of turbulent boundary layers that greatly impact and alter hydrodynamics in larger-scale CSAF units. CFD has been used to exploit this phenomenon to design a novel and scalable third-generation CSAF device that produces shear stresses well above the minimum required to inhibit cell deposition at a rotational frequency of only 20 rpm. This greatly reduced  allows scale-up to an effective filtration area of ca. 6 m2, which represents more than a 1,000-fold scale increase over the original CSAF design and would permit complete processing of 1,000 L of culture fluid in 2 hours.  64  Chapter 4  4 Characterizing solute dispersion in stacked membrane chromatography under flow-through conditions using a zonal rate model 4.1  Synopsis  Having developed a computational description of the hydrodynamics of the CSAF rotor chamber a characterization of the solute mass transport and binding kinetics within the associated membrane chromatography column is required for a comprehensive understanding of the physical phenomena within the entire technology. As discussed in Section 1.5 a necessary component of establishing a model for stacked membrane chromatography is the accurate accounting of contributions to band spreading occurring in the pre- and post-column hold-up volumes. Here, a new model for membrane chromatography that specifically accounts for non-uniform radial flow distributions within the inlet and outlet manifolds is derived and evaluated. This general model provides inherent flexibility, as it can be applied both to the stand–alone stackedmembrane chromatography capsule of the connected CSAF configuration (Figure 1.7) as well as the membrane stack associated with the original CSAF configuration (Figure 1.6). The model implementation and results reported in this chapter and in Chapter 5 are tailored to the connected configuration, as the resulting model equations are directly applicable to both CSAF technology and conventional membrane chromatography systems. The zonal rate model (ZRM) radially partitions the external hold-up volumes and membrane stack into virtual zones.  This partitioning allows for the independent  treatment of solute dispersion within the hold-up volumes of each zone and therefore  65  offers the potential to better capture the effect of the variable path lengths illustrated in Figure 1.9. The ZRM is first compared to data to determine the zonal configuration that best describes the dispersion of a model protein (ovalbumin) through the commercially available Mustang Q XT5 stacked-membrane anion-exchange chromatography module. The model is then used to decompose the observed dispersion to determine the magnitude of each contribution. All experiments and associated model calculations were performed under flow-through conditions to eliminate dispersion effects related to protein binding and desorption that might confound efforts to characterize and model non-idealities associated with the column hold-up volumes. These results therefore build on recent work performed on a simpler flat-sheet membrane unit for which an analytical solution to the tank network was possible, in part because dispersion in the membrane was neglected [166].  4.2  Theory  4.2.1 ZRM model description In its most general form, the ZRM radially partitions a stacked-membrane chromatography column (hold-up volumes and membrane stack inclusive) into any number n of zones as illustrated in Figure 1.8. There are then n membrane zones preceded by n feed distributor zones and followed by n eluent collector zones. The membrane stack is considered homogeneous in its properties, but is subject to different initial conditions at the membrane inlet within each zone due to differences in solute residence times within the associated zones of the feed distributor. Mobile and stationary phase composition therefore depend on the zone in addition to their typical dependence on axial position and time. Configuration of the zones can be achieved in a number of ways. For the illustrative case where n is taken to be 3, Figure 4.1a shows one possible arrangement of the virtual holdup zones for modules that possess both axial and radial flow components within the inlet and outlet manifolds. Radial flow in this configuration, hereafter referred to as the "standard" configuration, is modeled by forcing solute molecules in the feed that  66  Membrane Stack  Distributor  3a  PFR  Q  2a  1a  r uto  b i 3a str  B i t/ D e l n  I  Q  in  2a  Q3  Q2  Q1  Figure 4.1  Zone 2  Zone 1  Q2  Q1  3b 2b  1b  Membrane Stack Q3  Q2  PFR 1a  Zone 3  Q3  Q1  Zone 3  Zone 2  Zone 1  Q3  Q2  Q1  Collector  A  3b 2b  Co lle cto r  /O  out  utl e  t  1b  Stirred tank and membrane network representation of a stacked-membrane chromatography module with (A) standard and (B) diagonal flow distribution and collection.  67  physically travel through the outer radial region of the membrane to virtually pass through a common PFR and then inlet CSTRs 1a, 2a, and 3a. The total average residence time of a solute molecule traversing this path is therefore modeled as the sum of the average residence times for the PFR and the three inlet CSTRs, while that for a solute molecule traveling down the column centerline is given by the sum of the residence times for the PFR and CSTR 1a only. Flow paths and average residence times in the eluent collection manifold are modeled in an equivalent manner. Different configurations of this three zone model are possible, as illustrated by the “diagonal” configuration shown in Figure 4.1b, where solute molecules traveling through the outer radial region of the membrane now pass through a common inlet zone and then zone 3a. This configuration has the advantage that it provides a solute flow representation that more closely resembles the true flow path by explicitly accounting for sources of solute dispersion upstream of the feed distribution manifold and providing a unique hold-up volume for each zone. However, it requires an additional parameter to do so, which unnecessarily complicates the model, and therefore the standard ZRM configuration in Figure 4.1a was first focused upon with the aim of minimizing both the number of regressed parameters and overall model complexity. The required number of zones is optimized for the membrane chromatography module under investigation and should reflect a good compromise between complexity and accuracy of the model. An insufficient number of zones will violate the assumption that concentrations are radially homogeneous within a given zone, whereas an excessive number of zones results in insignificant concentration differences between the zones and an unnecessary complexity to the model. Each virtual hold-up zone is modeled as a CSTR using the standard mass balance relation: dc CSTR dt    c  1  CSTR in   c CSTR    (4-1)  68  where  is the residence time of the tank and is equal to the ratio of the tank volume VCSTR to the volumetric flow rate through the tank Qi ( = VCSTR /Qi); cinCSTR and c CSTR are the solute concentrations in the tank inlet and outlet, respectively. A PFR modeled as:  0 PFR t    PFR cout cin  t  t lag t  t lag  (4-2)  is placed in series within the flow network and accounts for an additional time lag that is uncoupled from the system dispersion. Here, the time lag, tlag, is equal to VPFR/Q. It should be noted that since the PFR provides a time lag but no dispersion, its relative position in the flow network is not important. Moreover, two or more such dead volumes within the network can be modeled with a single PFR. In the ZRM, the total time lag in the system is therefore accounted for by a single PFR placed before the tank/membrane network. The membrane zones are modeled as isolated membrane stacks of defined cross-sectional area embedded in the network of stirred tanks. All experiments and model calculations reported in this chapter are for non-binding conditions, and solute convection and axial dispersion within a membrane stack are therefore modeled in the traditional manner with a simplified form of the one-dimensional continuity equation for chromatography (Equation 1-13) where protein adsorption is neglected [66,70,149]:  c c  2c  u int  Dx 2 t z z  (4-3)  Here, c(z,t) is the solute concentration in the mobile phase at time t and distance z from the column inlet, Dx is the axial dispersion coefficient, and uint is the interstitial velocity given by the ratio of the superficial velocity usup to the membrane porosity m. Danckwerts’ boundary conditions are applied at the inlet (z = 0) and outlet (z = L) of each membrane stack [167]: u int cin t   u int c0, t   Dx  c z  (4-4) z 0  69  c z  0  (4-5)  z L  Initial conditions assume a step concentration increase from cin = 0 to cin = co at the root inlet of the tank network at time t = 0. As Equations 4-1 to 4-5 apply to flow-through conditions, they can be applied to any solute component introduced into the network. The flexibility inherent in partitioning both the external hold up volumes and column volume itself allows the ZRM to be applied to a variety of flow phenomena, including those found in commercial membrane chromatography modules and the prototype CSAF technology. In the membrane chromatography system modeled here, each zone within the membrane stack is treated identically and all flow non-uniformities are the result of the external volumes, but this does not necessarily have to be the case. Hydrodynamic, structural and mass-transfer properties within the column can be treated as zone dependent. For example, axial velocities and/or porosities can differ among column zones to account for the type of radial heterogeneity of interstitial velocities sometimes seen in packed bed chromatography due to wall effects [168-170]. 4.2.2 Computational methods The complete zonally segregated upstream hold-up volume, membrane stack and downstream hold-up volume network could be solved in a standard sequential manner by first solving the upstream hold-up volumes, then solving the membrane zones one by one with the upstream results giving the required initial conditions, and so on. However, that approach requires adaptive time stepping of the differential equation solver that necessarily complicates treatment of the initial and boundary conditions.  The entire  system of coupled model equations is therefore solved simultaneously in a novel and much more computationally efficient manner as follows. In the forward problem, where a breakthrough curve is computed for a given set of model parameters and specified column operating conditions, the partial differential equations for the membrane zones are first discretized along the axial coordinate by the method of lines [171]. The membrane zones are then coupled with the hold-up volume network as  70  described in Appendix D, resulting in a large system of differential equations (ODE). The MATLAB solver ode15s [172] is used for time integration and must be restarted at any step discontinuity in the boundary conditions of the root tank. With default settings the ode15s solver is predominantly occupied with numerically computing the system Jacobian, which is a matrix of size m by m that contains the partial derivatives of the right hand side of the nonlinear ODE system with respect to the m state variables of the spatially discretized model. This matrix is computed several orders of magnitude faster and more accurately by explicit differentiation of the corresponding model equations. Performance of the ode15s solver is further improved by a variable transformation that yields a diagonal mass matrix on the left hand side of the ODE system. In the inverse problem, unknown model parameters are estimated from chromatogram data by repeatedly solving and updating the entire equation system using MATLAB’s iterative optimization function lsqnonlin that uses an algorithm based on the interiorreflective Newton method [173,174]. Local minima are avoided by a multi-start strategy with varied initial values.  Unknown parameters in the ZRM that must be regressed  include those appearing in the PFR model (tlag) and possibly the continuity equation (Dx) if a value is not already available. The ZRM as configured in Figure 4.1a also requires values for two residence times per zone. This number can be halved if the chromatography module can be treated as axially symmetric, meaning the feeddistribution tank and eluent-collection tank within a given zone have identical average residence times (i.e. 1a = 1b, 2a = 2b, etc.). Additionally, as derived in Appendix D, a flow fraction i value is required for each additional zone beyond zone 1. The need to develop the highly efficient solution algorithm described above arises because a typical system with one solute component, three zones, and 1000 axial knots is described by a system of m = 1·2·3·1000 = 6000 equations, and one with three components, four zones and 1000 axial knots yields m = 3·3·4·1000 = 36,000 equations. Since the computational effort increases nonlinearly with both the number of equations and unknown model parameters, a high computation speed is crucial for the practical  71  applicability of the ZRM. The implementation of the computational algorithm described above enables a 2.20 GHz desktop computer to solve the forward problem in ca. 10 s and the inverse optimization of 4 unknown parameters in less than 5 min.  The set of  MATLAB programs that together constitute the ZRM are shown in Appendix E. Moment analysis was performed numerically on the model simulated breakthrough curves using MATLAB’s built-in numerical integrator trapz which uses the trapezoidal method. The programming code used in this analysis is found in Appendix F.  4.3  Experimental materials and methods  Mustang Q XT5 anion-exchange membrane chromatography modules (generously donated by Pall Inc. (East Hills, NY, USA)) containing modified hydrophilic polyethersulfone membranes were used in this study. The height and cross-sectional area of the stacked membrane bed in these modules are 2.20 mm and 22.06 cm 2, respectively, and the hold-up volume is 3.21 mL for both the feed-distribution and eluent-collection manifolds. The nominal pore size and porosity m of the membranes is 0.8 m and 0.70 ± 0.05 based on data provided by the manufacturer. The XT5 capsule was attached to an AKTAexplorer FPLC system controlled by Unicorn 4.12 software (Amersham Pharmacia, Uppsala, Sweden). Breakthrough experiments were performed under non-binding conditions using 1 g L-1 ovalbumin (Sigma, Oakville, Canada) and a running buffer of 1 M NaCl in 25 mM Tris, pH 8.0 (Fisher Scientific, Ottawa, Canada). Experiments were performed at flow rates of 1.5, 5, 10 and 20 mL min-1. The membrane stack was cleaned following the suggested protocol of the manufacturer using a 1 N NaOH solution followed by a wash with 1 M NaCl in 25 mM H3PO4 (Fisher Scientific, Ottawa, Canada).  4.4  Results and discussions  4.4.1 Comparison of existing hold-up volume models with experiment Breakthrough curves under non-binding conditions for frontal loading of ovalbumin at four distinct flow rates (1.5, 5, 10 and 20 mL min-1) were measured and first used to 72  evaluate the performance of both a “base-case” stacked-membrane chromatography model that ignores all hold-up volume effects, as well as the more advanced “one-zone” model proposed by Roper and Lightfoot [136] that includes terms for estimating solute dispersion in the feed distributor (linear PFR to CSTR sequence), membrane stack (1D rate theory of chromatography), and eluent-collection manifold (linear CSTR to PFR sequence). Figure 4.2a compares base-case model results to breakthrough data at 1.5 mL min-1 for the case where Dx is set to 7 x 10-11 m2 s-1, the diffusivity of ovalbumin in water at the column operating temperature [175]. This model calculation follows from the fact that the axial Peclet number Peax is generally quite large in stacked-membrane chromatography, and several studies have shown for such systems that the axial diffusion coefficient can be used in place of the axial dispersion coefficient with minimal error [e.g. 70]. Thus, Dx is expected to be near 7.0 x 10-11 m2 s-1, and this value is likewise in good agreement with available experimental values for membrane and monolithic columns (where, typically, 10-12 ≤ Dx ≤ 10-9 m2 s-1) [133,135] and appropriate correlations for estimating Dx from the physical structure and operating conditions of the column [176]. As shown in Figure 4.2a, the base-case model shows essentially no agreement with the breakthrough data when Dx is set equal to 7.0 x 10-11 m2 s-1, predicting almost no band broadening. As the experiments were conducted under non-binding conditions and low solute (ovalbumin) concentrations, the significant discrepancy between the base-case model and experiment cannot be attributed to dispersion effects related to slow kinetics of protein binding to or desorption from the stationary phase. The results therefore provide compelling evidence that solute dispersion within the inlet and outlet hold-up volumes makes a significant contribution to broadening of the breakthrough curve. Some amelioration of the base-case model error can of course be achieved by treating Dx as an apparent or lumped parameter that includes the contributions of the hold-up volumes to solute dispersion. This adaptation of the base case was therefore tested by globally regressing Dx to the complete set of breakthrough data. An improvement in model performance is observed (Figure 4.2a), but significant discrepancies remain. Moreover, minimization of model error results in a regressed value of Dx (1.0 x 10-4 m2 s-1) that is clearly unrealistic from the perspective of specifically treating axial dispersion within the membrane stack, as it is at least 5 orders of magnitude higher than the 73  Figure 4.2  Model fit of experimental breakthrough curves (●) for loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 module under non-binding conditions at a flow rate of 1.5 mL min-1 using: (A) the base-case membrane chromatography model which is comprised of a plug flow reactor in series with the 1D rate theory of chromatography; the grey line represents base-case model results for a Dx of 7 x 10-11 m2 s-1 and the black line for a regressed Dx of 1 x 10-4 m2 s-1 (B) the one-zone membrane chromatography model comprised of a plug flow reactor, a CSTR, the membrane stack and a second CSTR in series; the grey line represents model results for a Dx of 7 x 10-11 m2 s-1 and the black line for a regressed Dx of 1 x 10-5 m2 s-1. 74  expected value of Dx for this system. The results in Figure 4.2b show that the one-zone model of Roper and Lightfoot substantially addresses many of the shortcomings of the base-case model. When the flow rate is 1.5 mL min-1 and Dx is set to 7.0 x 10-11 m2 s-1, the one-zone model clearly captures the early stages of breakthrough much better than the base-case model.  However, it does not accurately reproduce the later stages of  breakthrough, including the slow approach to saturation of the membrane interstitial volume and associated hold-up volumes. These model errors, which also include the prediction of an initial breakthrough trend that is too shallow when compared to experiment, are similar to what has been previously reported for other one-zone type models [177]. The same errors are also observed at higher flow rates, though their magnitude decreases with increasing flow rate as indicated by the sum of the square of residuals data reported in Table 4.1. A small reduction in model error can again be realized by treating Dx as a lumped adjustable parameter and regressing it to the complete set of breakthrough data (Figure 4.2b). As before, however, the "optimal" value of the regressed Dx (1.0 x 10-5 m2 s-1) is completely unrealistic as an indicator of solute dispersion specifically within the membrane stack.  Table 4.1  ∞  Sum of squared residuals for the fitting of various models to experimental breakthrough data at four different flow rates for loading of ovalbumin (1 g L-1) on a Mustang XT5 column under non-binding conditions. Dx was 7 x 10-11 m2 s-1 for all model calculations.  Flow rate (mL min-1)  Base case  1 Zone∞  2 Zones (w/ & w/o symmetry)  3 Zones (w/ symmetry)  1.5  9.1374  0.5173  0.0034  0.0029  5  12.7519  0.3198  0.0077  0.0077  10  13.7395  0.0302  0.0074  0.0058  20  14.8522  0.0212  0.0133  0.0079  Results computed using the Roper and Lightfoot model [136] with Dx = 7 x 10-11 m2 s-1  75  4.4.2 ZRM model configuration optimization and application Flow-through conditions are best suited for determining the structure of the hold-up volume network in the ZRM, as the unknown system parameters are then limited to the residence times for the zonal CSTRs, the lag time for the required PFR, and Dx if it is not assumed equal to the diffusivity of ovalbumin. Dispersion effects associated with protein binding and elution are avoided and all other required model parameters are known, including the stacked membrane thickness and porosity, as well as the hold-up volumes of the inlet and outlet manifolds. In the XT5 module, the geometries of those two manifolds are simple mirror images, and the manifolds therefore share the same hold-up volume, 3.21 mL each. As a result, the average solute residence time is expected to be the same for the two manifolds. A number of configurations of the ZRM was globally fit to the complete set of ovalbumin breakthrough data collected from the XT5 membrane chromatography module and then the sum of the square of residuals was compared. Both the number of zones and the validity of the manifold symmetry approximation were tested to determine the ZRM configuration that best describes the experimental data with the minimum number of regressed parameters. Table 4.1 reports sum of the square of residuals data for the fitting of several configurations of the standard ZRM to the breakthrough curves measured for frontal loading of ovalbumin at four different flow rates. All configurations of the ZRM tested resulted in significant improvement in model accuracy when compared to either the basecase or single-zone model. Indeed, as illustrated in Figures 4.3a and 4.3b, calculated breakthrough curves from the simplest ZRM configuration tested (two zones) essentially superimpose on the experimental data, and this result was true at all flow rates as evidenced by the sum of the squared residual values reported in Table 4.1. No change in ZRM performance was observed when the symmetry approximation was invoked, as should be expected from the geometry of the XT5 capsule. As a result, though the breakthrough data are also perfectly described by the ZRM configured either in two zones without symmetry or in three zones with tank symmetry (Table 4.1), the added  76  Figure 4.3  Symmetric two-zone ZRM fit (solid line) of the experimental breakthrough data (●) for ovalbumin (co = 1 g L-1) loaded on a Mustang XT5 module under non-binding conditions at a flow rate of (A) 1.5 mL min-1 and (B) 5 mL min-1. Dx was 7 x 10-11 m2 s-1 in all calculations.  77  complexity of those configurations is not required. The same arguments also made it unnecessary to evaluate the performance of the “diagonal” configuration of the ZRM shown in Figure 4.1b. Within the ZRM, the flow fraction through the outer zone, 2is independent of flow rate and was therefore directly regressed to the entire data set, giving a value of 0.26. This relatively low 2indicates that the ZRM can successfully address the deficiencies in the single-zone model to quantitatively capture solute dispersion in a stacked-membrane chromatography module by directing a relatively small percentage of the feed through zone 2. The values of the remaining parameters of the symmetric two-zone ZRM are reported in Table 4.2. All of these regressed times necessarily depend on flow rate since the associated hold-up volumes are constant. The product of the lag time tlag and Q is a constant characteristic dead volume (ca. 0.38 mL) in accordance with the solution of the PFR model (Equation 4-2). Similarly, the residence time for each CSTR has a linear dependence on 1/Q, as is required by the solution to Equation 4-1. From a modeling perspective, what is important is that the values of and the disparity between the inner and outer zone residence times decrease with increasing Q. Solute dispersion effects associated with unequal flow paths through different regions of the inlet and outlet manifolds are therefore predicted to be largest at low flow rates. All CSTR residence  Table 4.2  Parameters for the symmetric two-zone ZRM configuration regressed from experimental breakthrough data at four different flow rates for ovalbumin (co = 1 g L-1) loaded on a Mustang XT5 column under nonbinding conditions. Dx was set to 7 x 10-11 m2 s-1. Flow rate (mL min-1)  tlag (s)  inner  outer  (s)  (s)  1.5  15.25  45.4  125.3  5  4.57  20.5  37.6  10  2.29  11.5  18.8  20  1.14  6.1  9.4  78  times approach zero near the highest flow rates recommended for operation of the module. As a result, differences between the symmetric two-zone ZRM and the singlezone model are diminished (Table 4.1), indicating that the ZRM effectively reduces to the simpler one-zone model of Roper and Lightfoot (but not the base-case model) at high Q. 4.4.3 Decomposition of system dispersion using the ZRM In addition to accurate prediction of breakthrough, the ZRM provides the ability to decompose solute dispersion in order to analyze contributions taking place in each radial zone as a function of time and flow rate.  Figure 4.4 reports contributions to solute  breakthrough from the inner and outer radial zones computed by the symmetric two-zone configuration of the ZRM for two different feed flow rates. The trends are qualitatively similar between the two flow rates. Significantly greater solute dispersion is predicted within zone 2 due to the longer average residence time of solute passing through the outer radial positions of each manifold. Breakthrough of solute from zone 2 is therefore delayed relative to that from zone 1 and, as a result, zone 2 does not contribute to the initial breakthrough of solute from the column. The ZRM therefore predicts a relatively steep initial breakthrough, as observed experimentally but not well described by previous models. The ZRM also predicts that when breakthrough from zone 1 is complete, the loading of zone 2 is less than half complete and the remaining slow approach to system saturation is fully determined by dispersion processes within zone 2 during the second half of its loading. For the 1.5 mL min-1 flow rate, the variance (v) of the various dispersion processes within each zone are reported in Table 4.3 along with the total variance generated by the zone. For either zone, solute dispersion within the feed-distribution and eluent-collection manifolds accounts for at least 98% of the total dispersion. As should be expected for the commercial ion-exchange membranes used in this study, the membrane stack is highly efficient, offering a height equivalent to a theoretical plate, HETP, of between 1 and 4 µm over the range of flow rates tested and contributing very little to the observed solute dispersion within the module. Indeed, the near step-change breakthrough curve (grey  79  Figure 4.4  Contributions of the inner and outer radial zones to elution band spreading of ovalbumin (co = 1 g L-1) loaded on the XT5 module at (A) 1.5 mL min-1 or (B) 20 mL min-1. The symmetric two-zone configuration of the ZRM was used for all calculations. The residence times of the inner and outer tanks were 45.4 and 125.3 s, respectively, and the time lag was 15.25 s. Dx was set to 7 x 10-11 m2 s-1.  80  Table 4.3  Breakdown of ZRM predicted solute elution variance v (s) for the 1.5 mL min-1 breakthrough data as a function of flow zone and flow elements within the zone. Source of Variance  Zone 1  Zone 2  Total  81.45 s  267.9 s  Feed distribution manifold  48.68 s  141.6 s  Membrane stack  1.62 s  0.33 s  Eluent collection manifold  31.15 s  126.0 s  line) in Figure 4.2a computed using the base-case model represents the degree of elution band broadening that would be expected if the external hold-up volumes provided no contribution to solute dispersion. It is interesting to note that the contribution of the membrane stack to the total solute dispersion is smaller in zone 2 (v = 0.33 s) than in zone 1 (v = 1.62 s) despite the fact that both Dx and uint are the same within the two zones. This is due to differences in the shape and slope of the concentration gradient of solute entering the membrane stack from the two different zones of the feed distribution manifold. Due to the higher dispersion of solute exiting zone 2 of that manifold, d2c/dz2 in Equation 4-3 takes on smaller values in general, resulting in lower dispersion within zone 2 of the membrane stack. Differences in initial conditions likewise account for the unique variance contributions made by each CSTR within the network. The ZRM can also be used to explore operational space to identify system configurations or conditions where improvements in system performance might be realized.  For  example, Figure 4.5 shows ZRM predictions of the effect on total solute dispersion, quantified here by the computed normalized second central moment (  2 / 1  ) where 2   2 equals v2 and 1 is the first absolute moment, of changes to the values of 2 and the column Peclet number Peax. When Peax > 100, which is true at 1.5 mL min-1 (Peax = 510) and typically the case in membrane chromatography, the contribution of the membrane to  81  Figure 4.5  Normalized second moment of ZRM simulated breakthrough curves as a function of axial Peclet number Peax and flow fraction through the outer zone (2) for the case where the column length is 2.2 x 10-3 m and the flow rate is (A) 1.5 mL min-1 or (B) 20 mL min-1. The symmetric twozone configuration of the ZRM was used for all calculations.  82  solute dispersion is negligible and the ZRM predicts that a sharpening of solute breakthrough can be achieved by engineering a set of inlet and outlet manifolds that serve to reduce 2. As shown in Figure 4.5, even a small decrease in 2 could provide a significant improvement, as the ZRM predicts that  2 / 1  decreases rapidly with 2  decreasing 2 when Peax > 100. The results in Figure 4.5 also show that solute dispersion within the membrane stack is predicted to contribute significantly to broadening of the elution band when Peax drops to a value near 10 or lower. Thus, at a flow rate of 1.5 mL min-1, a 2.2 mm thick membrane stack would contribute to the total variance if it were characterized by a Dx of 10-9 m2 s-1, which is at the high end of experimental Dx values reported for membrane and monolithic columns [133,135,178]. In this regime (Dx ≥ 10-9 m2 s-1), band broadening is predicted to depend strongly on Dx but, because the contributions of the inlet and outlet manifolds are reduced, show a much weaker dependence on 2. Operationally significant dispersion within the column can also be realized through an increase in column length. As shown in Figure 4.6, an order of magnitude increase in the membrane bed height to 2.2 cm is predicted to increase  2 3-fold (note the log scale), with further increases in bed height resulting in more significant increases in band broadening due to the membrane stack then making the dominant contribution to v. When configured properly, the ZRM can therefore provide guidance on effective ways to increase column capacity by predicting the effects of increasing both bed diameter (which may increase hold-up volumes, associate residence times and number of zones) and bed height on column loading efficiency. Finally, as shown in Section 4.4.2, the symmetric two-zone ZRM reduces to the one-zone model at high feed flow rates. In terms of the ZRM, this would mean that the total variance loses its dependence on 2 with increasing flow rate. This is indeed the case (Figure 4.5b). At a feed flow rate of 20 mL min-1,  2 / 1  becomes nearly insensitive 2  to 2.  83  Figure 4.6  4.5  Second central moment of ZRM simulated breakthrough curves as a function of column length and flow fraction through the outer zone (2) for the case where the flow rate is 1.5 mL min-1 and Dx is 7 x 10-11 m2 s-1.  Concluding remarks  In stacked-membrane chromatography the large column diameter to length ratio naturally leads to variations in solute residence times within the feed-distribution and eluentcollection manifolds.  The ZRM accounts for these unequal flow paths by radially  partitioning the membrane stack and associated hold-up volumes into virtual zones. By fitting the model to a small set of breakthrough data obtained under non-binding conditions - thus avoiding any dispersion effects associated with protein adsorption and desorption - an appropriate ZRM structure can be identified and then used both to quantitatively match breakthrough data over a range of operating conditions, including those flow rates where previous membrane-chromatography models fail, and to decompose the breakthrough data to identify and possibly correct major sources of band broadening within the system. The results highlight the importance of proper design of the external hold-up volumes, especially in systems where the membrane volume is equal to or less than the external hold-up volumes. 84  In the following chapter, the ZRM is used to characterize stacked-membrane chromatography under binding conditions to permit the regression of binding kinetics parameters from limited data, the prediction of breakthrough curves, and the optimization and scale-up of a separation from data acquired on the scaled-down XT5 module. Treating the hydrodynamic and adsorption effects separately will enable scale-up to proceed in an economical and efficient manner. The zonal arrangement as well as tank sizes and flow ratios must be defined for different scales and geometries following the procedure outlined in this chapter. However, the binding parameters determined from a scaled-down device are directly transferable across scales and geometries as they are independent of the flow effects in the external hold-up volumes. This greatly reduces the costs associated with optimizing large-scale chromatography modules.  85  Chapter 5  5 Characterizing stacked membrane chromatography under protein retention conditions using a zonal rate model 5.1  Synopsis  In the previous chapter the novel ZRM was developed to accurately account for solute dispersion within membrane chromatography modules, including pre- and post-column dispersion resulting from radial flows within the hold-up volumes.  The ZRM  quantitatively captures this non-ideality by expanding on sequential PFR/CSTR models [131,133,136,142,143,149,150] through the partitioning of the membrane stack and its associated pre- and post-column hold-up volumes into virtual zones that account for differences in radial path lengths traveled by solute molecules within the feed and eluent distributors of the column. A defined set of pulse or frontal loading experiments under non-retention conditions were used to determine the extra-column parameters and zone configuration of the ZRM.  Through careful accounting of contributions to band  spreading occurring in the pre- and post-column hold-up volumes the ZRM provides a framework to properly isolate and characterize mass transport and intrinsic binding within the membrane stack itself. To date, protein binding to membrane adsorbers has most often been modeled using standard Langmuir kinetics [66,70,133,135,142,143,149,179]. However, when combined with appropriate representations of the GRM, Langmuir theory cannot capture the complex elution behavior described in Section 1.5, and this has motivated a number of researchers to seek more realistic binding models that might serve to improve predictions and advance knowledge of protein retention mechanisms in IEX membrane (IEXM) chromatography systems. For example, Sarti and coworkers [131] greatly improved BTC modeling by implementing bi-Langmuir adsorption kinetics. Their approach therefore  86  assumes that heterogeneity in the binding energy of the sorbent accounts for the complex elution behavior. Alternative theories have been advanced by Yang and Etzel [143], who showed that BTC asymmetry in IEXM columns can also be described using either a steric-hindrance/SMA type model to account for the shielding of binding sites by large protein molecules or a spreading type model that considers the conformational changes a protein may undergo when bound to the sorbent. However, each adsorption model typically accounts for only one of the adsorption non-idealities noted in Section 1.5, in large part because such simplification is required to reduce the complexity and parameterization of the model so as to permit accurate fitting to a limited set of data, in this instance typically obtained on-column. Here the “non-retention” version of the ZRM (nr-ZRM) is extended to describe retention and breakthrough of ovalbumin frontally loaded on the Mustang Q XT5 chromatography module under binding conditions. Dimensionless groups are determined to reduce the GRM to its most appropriate form, and isotherm and breakthrough data are used to identify both an appropriate binding model and the dominant non-idealities associated with the binding process.  Appropriate forms of Langmuir, bi-Langmuir, SMA and  spreading theory are evaluated for their ability to fundamentally describe and quantitatively predict BTCs when implemented within the framework of the ZRM model. As part of this process, the contribution of film mass transport to protein breakthrough behavior is also investigated.  5.2  Theory  5.2.1  Zonal rate model  The basic structure of the ZRM is described fully for non-retention conditions in Section 4.2.1 and is shown in Figure 5.1 for the specific case where non-uniform flows and dispersion within the extra-column volumes are modeled with 2 radial zones. To extend the nr-ZRM to describe retention and breakthrough of protein under binding conditions Equation 4-3 is replaced with the one-dimensional continuity equation for a retained species [66,70,134]:  87  Figure 5.1  Stirred tank and membrane network representation of a stacked-membrane chromatography module with two radial zones. The total flow rate, Q, is split between zone 1, Q1, and zone 2, Q2. The stirred tank network is symmetric in that the residence times, i, of the two tanks in each radial zone are equivalent.  c c  2 c 1   m  s  u int  Dx 2  t z  m t z  (5-1)  which is solved for each zone in time and space as before with Equations 4-1 (or Equation D-1 in the case of a tank with multiple inlets, see Appendix D) and 4-2. In Equation 5-1, s(z,t) is the total solute concentration in the stationary phase of the zone. As in the nr-ZRM, Danckwerts’ boundary conditions (Equations 4-4 and 4-5) are applied at the column inlet (z = 0) and outlet (z = L) of each zone. [167]. In Equation 5-1, the volume of the stationary phase comprises both the volume of the solid sorbent and that of the liquid film that surrounds the sorbent. The sink term is therefore given by  c f s q 4 f  1   f   k f c  c f   k f c  c f t t t d p    (5-2)  where f is the fraction of the total stationary phase volume occupied by the liquid film, cf is the solute concentration within the film, q is the sorbate concentration, kf is the film  88  mass-transfer coefficient, and k f is (4/dp)kf where dp is the nominal pore diameter of the membrane and 4/dp gives the pore surface area to volume ratio. Equation 5-2 solves for the time derivative of the solute concentration in the film, cf, when the time derivative of q is given by the binding model. This careful treatment of film mass transfer was not considered in previous models of membrane chromatography [61,70,131,138,143,149,151,177]. However, it is both rigorously correct and very useful, in part because it properly partitions the column volume into three distinct volume fractions: the fractions occupied by the mobile phase m, the stagnant film (1-m)f, and the stationary phase (1-m)(1-f), respectively. Technically the column volume could be partitioned differently, but the volume fractions used here offer the advantage that the resulting model equations are analogous to those of the linear driving force model for packed-bed chromatography [179]. 5.2.2 Intrinsic binding rate models Solution of Equations 5-1 and 5-2 requires an appropriate model describing the intrinsic kinetics of protein binding to and desorption from the sorbent.  Models that have been  used to describe membrane adsorbers and will be examined here include the traditional Langmuir model [e.g. 61,70,135,142], the bi-Langmuir model [131], the SMA model [147], and a general form of the spreading model [143,145,146]. Protein binding rates are described in the traditional Langmuir model as:  q  k a c f qm  q   k d q t  (5-3)  where ka and kd are the adsorption and desorption rate constants and qm is the saturation capacity. This simplest of binding isotherm models assumes sorbate molecules are noninteracting and that all binding sites on the sorbent are identical and independent. In bi-Langmuir theory, where the sorbent is modeled as offering two energetically distinct types of independent binding sites, protein binding kinetics are described by:  89  q q1 q 2   t t t  (5-4)  q1  k a ,1c f q m,1  q1   k d ,1 q1 t  (5-5a)  q 2  k a , 2 c f q m, 2  q 2   k d , 2 q 2 t  (5-5b)  where q1 and q2 represent the concentrations of solute bound to sites of type 1 and type 2 respectively; ka,1, kd,1 and qm,1 are the binding parameters associated with the high-energy binding site; and ka,2, kd,2 and qm,2 are the binding parameters associated with the lowerenergy binding site. The SMA model accounts for both the effect of steric shielding and the effect of co-ion (A) concentration on the adsorption of a protein macro-ion (P) of characteristic charge v onto a sorbent (S). It is generally applied in its equilibrium isotherm form, but for protein binding reactions of the form:  P v  v A  S   vA  P  S  (5-6)  it also yields intrinsic binding rate equations of the form:  q v  k a c f q A  k d qc Av t  (5-7)  q A    (v   )q  (5-8)  where   is the ion capacity of the stationary phase, q A is the concentration of bound salt coions that are available for exchange with the protein and  is the steric factor. Finally, Clark and co-workers [145] have developed a general spreading model in which a protein molecule in the adsorbed layer can be reversibly transferred among a continuum of bound states, each characterized by a unique affinity and projected area of interaction 90  due to changes in the orientation or conformation of the bound protein [180-182]. Bound proteins at any state can also transfer to and from the bulk fluid. For the case where two adsorbed protein states (1 and 2) are considered, the general spreading model reduces to the following rate equations:  q q1 q 2   t t t  where  q1  (5-9)  q1  ka ,1c f  k12q1 qm,1  q1  q2   kd ,1q1  k21q2 t  (5-10a)  q2  ka , 2c f  k12q1 qm,1  q1  q2   k21  kd , 2 q2 t  (5-10b)  and  q2  represent  the  concentrations  of  bound  protein  in  orientational/conformational states 1 and 2 respectively,  is the ratio of the sorbent surface area occupied by a bound protein in state 2 relative to that in state 1, and the rate constants are as defined in Figure 5.2. In Equations 5-10, qm,1 represents the hypothetical saturation capacity of the stationary phase if all adsorbed proteins are in state 1. The true  Figure 5.2  Schematic of the mechanism of adsorption of a protein molecule according to the simplified spreading model, where protein molecules can reversibly adsorb to a homogeneous binding surface in two different conformations/orientations.  91  qm is given by q1 + q2 at saturation conditions. This model is similar to one used previously to model reorientation of adsorbed proteins on a solid surface [180], but takes into account the surface area of bound molecules as well as the dependence of spreading on surface coverage. It reduces to the spreading model developed by Lundstrom [146] and implemented in membrane chromatography by Yang and Etzel [143] when the spreading of the adsorbed protein is considered irreversible and protein from the bulk cannot adsorb directly to state 2 (i.e. when k21 = ka,2 = 0). Each of these models is based on a distinctly different mechanism of binding, yet all have been used to describe protein retention in membrane chromatography columns. The goal is to establish a framework that permits more rigorous evaluation of the merits of each model, particularly with respect to providing a realistic interpretation of the intrinsic binding process within the model system. Note that since the protein binding reaction in the membrane chromatography system is conducted under flow, qm and qm,i in the rate models derived above are typically replaced to improve model fit with their corresponding dynamic (*) saturation capacities, q m and q m ,i , respectively.  As Q  approaches zero, q m and q m ,i approach qm and qm,i, respectively. Above some critical flow rate, q m and q m ,i will decrease with increasing Q.  5.2.3 Computational methods A new highly efficient computation strategy was developed to solve the ZRM and is described in detail for the case where the protein is not retained in Section 4.2.2. Here, the added features of the computational algorithm when the ZRM is applied to protein loading are described. When the partial differential equations describing mass transport are coupled with the adsorption kinetic equations of the SMA isotherm a large system of differential-algebraic equations results. This system, similarly to systems of ODEs, is solved simultaneously to synchronize the adaptive time stepping of the differential equation solver.  The  MATLAB solver ode15s [172] is used for time integration and, in the case of algebraic  92  equations, this solver must be restarted at any step discontinuity in the boundary conditions of the root tank. In the inverse problem, as in the case of the nr-ZRM, unknown model parameters are estimated from chromatogram data by repeatedly solving and updating the entire equation system using MATLAB’s iterative optimization function lsqnonlin. Unknown binding parameters that need to be regressed from the data include the rate constants, the binding capacities, the axial dispersion coefficient and, in the case where film diffusion is non-negligible, the stationary phase ratio f and the film mass-transfer coefficient kf. Universal regressions were performed simultaneously across all flow rates as some parameters, such as the adsorption and desorption rate constants, are flow rate independent whereas others, such as the rate constants associated with transformation between bound states as well as film diffusion coefficients, vary with flow rate.  5.3  Experimental materials and methods  The experimental materials and equipment are described in Section 4.3. Breakthrough experiments were performed under binding conditions using 1 g L-1 ovalbumin (Sigma, Oakville, Canada). Frontal analysis at 0.5 mL min-1 for a ladder of ovalbumin feed concentrations ranging from co = 0.2 to 10 g L-1 was used to determine equilibrium isotherm and (static) saturation capacity data under column loading conditions. The equilibrium concentration of bound protein q was computed from the BTC data as   q   co  c Q Vcolumn dt , where Q and Vcolumn are the feed flow rate and column 0  volume, respectively. The load and elution buffers were 25 mM Tris-HCl (pH 8.0) and 1 M NaCl/25 mM Tris-HCl (pH 8.0) respectively (Fisher Scientific, Ottawa, Canada) and were filtered and degassed under vacuum prior to use. The membrane capsule was cleaned as per the manufacturer’s suggested protocol using a 1 N NaOH solution followed by a wash with 1 M NaCl in 25 mM H3PO4 (Fisher Scientific, Ottawa, Canada).  93  5.4  Results and discussions  5.4.1 Dimensionless group analysis In membrane chromatography, solute diffusion through the hydrodynamic boundary layer at the sorbent surface can contribute to the overall resistance to protein mass transfer [133,135,142,148]. It is, however, typically neglected in membrane chromatography models [61,70,131,138,143,149,151,177]. The justification for doing so is generally made by comparing the characteristic time for film mass transfer (tfilm) to that for axial convection (tconv) through the entire membrane stack [70,131,143,149,151,177]:  t film   d p2 4 Dx    L  t conv u int  (5-11)  However, as the intent of this analysis is to determine limiting rates to local equilibrium, it is better made by comparing tfilm to tplate, the characteristic time for solute convection through the interstitial volume of a single theoretical mass-transfer unit. This ratio follows from non-dimensionalizing Equations 5-1 and 5-2 using the height equivalent to a theoretical plate, HETP, as the characteristic length scale:  c~ c~ 1  2 c~ 1   m  1 ~ ~ c  c f    ~  t ~ z Pe ax ~  m Pe f z2    (5-12)  z HETP  (5-13)  The non-dimensionalized variables are then given by:  c=  c , co  tu int ~ t  HETP  and  z=  Similarly, the rate equations for intrinsic protein adsorption to the membrane can be nondimensionalized with respect to dp and kf.  Thus, the resultant expression for the  Langmuir model is:  kd d p ~ q~ 1 ~  c 1  q~   q tˆ Dk f kf  (5-14)  94  where the non-dimensionalized variables for the stationary phase are given by:  q=  Three  dimensionless  groups  q qm are  and  thus  tˆ =  tk f  (5-15)  dp  defined:  an  axial  Peclet  number,  Pe ax  uint HETP Dx , a film Peclet number, Pe f  uint d p 4HETPk f , which gives the  ratio of characteristic times for convective and film mass transfer, and a Damköhler number, Dk f = k f k a co d p , which gives the ratio of characteristic times for solute adsorption and film diffusion. All numbers represent appropriate metrics for determining the relative contribution of film mass transfer to band broadening. HETP values were determined (Figure 5.3) from experimental BTCs for frontal loading of retained ovalbumin at four distinct flow rates (1.5, 5, 10 and 20 mL min-1). The film mass transfer coefficient was estimated from the correlation of Athalye et al. [183]: 3  1.09  d p  m u int  Dm   64    kf  dp   m  Dm      1/ 3  (5-16)  which has previously been used in membrane chromatography [70]. Dm is the molecular diffusivity, 7 x 10-11 m2 s-1, as computed by the method of Tyn and Gusek [175]. This analysis permits Pef to be estimated as a function of linear velocity. As shown in Figure 5.3, over the flow conditions recommended for the capsule, Pef ranges from 3.2 x 10-4 to 1.5 x 10-3, indicating that the resistance to film mass transfer is negligible and can safely be neglected with no loss of model accuracy. This yields a simplified form of the model in which Equations 5-1 and 5-2 reduce to:  c c  2 c 1   m  q  u int   m Dx 2  t z  m t z  (5-17)  95  Figure 5.3  Calculated height equivalent to a theoretical plate (HETP), estimated film Peclet number (Pef) and calculated axial Peclet number (Peax) for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 column under binding conditions. HETP values are based on plate theory and derived from breakthrough curves at 1.5, 5, 10 and 20 mL min-1. Film diffusion coefficients are estimated from Equation 5-16.  and cf can be replaced with the concentration of solute in the bulk, c, in this and all intrinsic binding rate equations. Later it will be shown that use of Equation 5-17 is likewise supported by the value of Dkf. 5.4.2 Benchmark study – the Thomas model Transport of ovalbumin within a Mustang Q XT5 capsule is characterized by a Dx of 0.9 (± 0.3) x 10-8 m2 s-1, which is typical of axial dispersion coefficients previously reported for globular proteins in stacked-membrane chromatography modules [135].  Peax  therefore lies between ca. 4 and 60 for the range of feed flow rates applied to the capsule (Figure 5.3). Based on theoretical considerations, Suen and Etzel [70] have argued that breakthrough from a stacked membrane column operated under loading conditions should obey the Thomas model when Peax ≥ 40 and film mass-transfer resistances can be 96  neglected (i.e. Pef << 1). The Thomas model, an analytical solution for Equation 5-17 when Dx = 0 and Langmuir kinetics are applied, is given by: J n r ,nT  c  c0 J n r ,nT   1  J n r ,nT exp 1  1 r n  nT   n  1   m qm k a L  T  (5-19)   m u int r = 1+  (5-18)  co Kd  (5-20)   m K d r  u int t   1  1   m qm  L   (5-21)  where the J function in Equation 5-18 is defined as:      J x, y   1  exp  y  exp   I o 2 y d x  0  (5-22)  and Io is the modified Bessel function of zeroth order. As it considers only convective mass transport and Langmuir binding kinetics, the Thomas model effectively represents the best loading and breakthrough performance that can be achieved in membrane chromatography. It therefore provides a useful metric against which to compare actual column performance. Figure 5.4 compares experimental BTCs for ovalbumin loaded at a feed concentration (co = 1 g L-1) lying within the linear region of the isotherm to those calculated by a benchmark model comprised of a linear sequence of a PFR, a pre-column CSTR, the Thomas model, and a post-column CSTR of equal volume to the first CSTR. The structural parameters for the PFR and CSTR elements were determined under nonbinding conditions and are tabulated in Table 5.1. The Langmuir rate parameters (in parentheses in Table 5.2) were regressed to the BTC data shown for a feed flow rate of 20 mL min-1, where Peax is above 40. As is typically observed in both traditional and membrane chromatography of proteins,  qm* depends on uint and is therefore  independently fit to each experimental BTC.  97  Figure 5.4  Comparison of the analytical solution of the Thomas model coupled with a pre- and post-column CSTR to experimental breakthrough data (♦) for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 module under binding conditions at flow rates of (A) 1.5 mL min-1and (B) 20 mL min-1.  98  Table 5.1  Structural parameters for (i) the symmetric one-zone model proposed by Roper and Lightfoot [136] comprised of a PFR in series with the membrane stack sandwiched between two CSTRs of equal residence time and (ii) the symmetric two-zone nr-ZRM. Model parameters were determined by fitting experimental data for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 column under non-binding conditions at flow rates of 1.5, 5, 10 and 20 mL min-1. In both models, all dispersion occurring in the dead volume modeled by the PFR is accounted for by increasing the solute residence time of the first CSTR in the network by a flow rate independent amount ′, so that 1/1 = 1/′ + Q1/VCSTR.  ∞  Table 5.2  Parameter  One-zone model∞  Two-zone nr-ZRM  VPFR  0.38 mL  0.38 mL  ′1  0.27 min-1  0.64 min-1  V1  3.11 mL  2.19 mL  V2  -  0.82 mL  2  -  0.26  Results computed using the Roper and Lightfoot model [136].  Regressed parameter values for the Langmuir adsorption model determined from fitting BTC data for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 column under retained conditions at four different flow rates (1.5, 5, 10 and 20 mL min-1). The Langmuir adsorption model was used within the symmetric two-zone ZRM. The rate constants in parentheses represent the Langmuir kinetics parameters determined from fitting the corresponding Thomas model to the BTC data at 20 mL min-1. Flow rate (mL min-1)  ka (L g-1 s-1)  kd (s-1)  1.5 5 10 20  qm (g L-1) 350  5.12x10-3 (5.51x10-3)  1.11x10-3 (1.11x10-3)  373 411 402  99  The poor agreement of the Thomas model with experiment at 1.5 mL min-1 is expected and reflects, in part, the importance of axial dispersion to band broadening under slower load conditions where Peax ~ 1. At 20 mL min-1, the Thomas model captures the BTC data well, though it shows a slightly steeper approach to saturation than observed experimentally. As the rapid mass-transfer conditions (Peax ≥ 40; Pef << 1) specified by Suen and Etzel are met at 20 mL min-1, this small discrepancy with experiment is most likely related to non-ideal processes occurring in the Mustang Q XT5 capsule that are not associated with solution-phase mass transfer within the membrane stack.  Those  processes may be separated into: i) flow non-idealities within the external hold-up volumes of the capsule; and, ii) features of the intrinsic binding reaction not accounted for in Langmuir theory.  The contributions of these possible non-idealities can be  independently analyzed and modeled by noting that non-idealities associated with protein binding are not present under flow-through conditions. 5.4.3 Accounting for external flow non-idealities using the ZRM For feed flows of 1.5 mL min-1 and 20 mL min-1, Figure 5.5 compares BTCs for ovalbumin measured under flow-through (non-binding) conditions with those predicted by Roper and Lightfoot’s [136] model comprised of a PFR in series with a membrane stack sandwiched between pre- and post-column CSTRs of equal residence time (i.e. Equations 4-1 and 4-2). Results from a symmetric two-zone version of the ZRM in which the residence times of the two tanks within a given zone are equal are also reported. Parameters for the two models are reported in Table 5.1.  Note that a small  amount of solute dispersion occurs in the extra-column dead volume modeled by the PFR. In either model, this dispersion is accounted for by increasing the solute residence time 1in the first CSTR in the network by a flow rate independent amount ′. Thus, 1/1 = 1/′ + Q1/VCSTR. The results show that both models work well at high feed flow rates. This is expected since the nr-ZRM collapses to the Roper and Lightfoot model at high feed flows due to the very short solute residence times within extra-column volumes.  The results also  100  Figure 5.5  Model fits of experimental breakthrough data (♦) for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 module under non-binding conditions at flow rates of (A) 1.5 mL min-1 and (B) 20 mL min-1using: (▬) the symmetric one-zone model comprised of a PFR in series with the membrane stack sandwiched between two CSTRs of equal residence time; (▬) the symmetric two-zone nr-ZRM.  101  confirm previous reports [133,142,149] that the use of a linear one-zone PFR/CSTR sequence to model solute dispersion within the external hold-up volumes is not exact at lower flows, with model errors increasing with decreasing flow rate. In contrast, the twozone ZRM accurately captures inhomogeneous flow effects in the external hold-up volumes of the Mustang Q XT5 IEXM column at all feed flow rates. As a result, external and internal solute dispersion are properly decoupled at all Q by the ZRM to provide a rigorous framework for analyzing the fitness of various binding models to describe protein loading and breakthrough. 5.4.4 Evaluating putative binding models 5.4.4.1 Langmuir model Figure 5.6 compares BTCs calculated by a symmetric two-zone ZRM utilizing Langmuir binding kinetics (Equation 5-3) to experimental data for frontal loading of ovalbumin at flow rates ranging from 1.5 to 20 mL min-1. Note that this calculation differs from the benchmark (Thomas model) calculations by the fact that both axial dispersion and multiple flow paths are considered.  As has been observed in other membrane  chromatography processes [149,154], particularly those involving ovalbumin [153], q m* increases with increasing flow rate within the low-flow operating regime. In accordance with most chromatographic processes, it then decreases with increasing flow when Q exceeds a critical value, in this case ca. 10 mL min-1. The best fit of the two-zone ZRM utilizing Langmuir binding kinetics captures this trend (Table 5.2). However, Langmuir theory treats every binding site as independent and equivalent. Thus, from a fundamental perspective, Langmuir theory cannot predict an increase in loading capacity with Q. It merely correlates to this trend when the regressed value of q m* is allowed to depend on Q. The ZRM with Langmuir kinetics is far more accurate than the benchmark Thomas model, in large part due to the accounting of axial dispersion at low flows, but also through the proper accounting of extra-column flow non-idealities.  Nevertheless,  agreement with experiment remains relatively poor (Table 5.3), indicating that the Langmuir binding model oversimplifies ovalbumin adsorption to the IEX membranes of  102  Figure 5.6  Comparison of the symmetric two-zone ZRM incorporating the Langmuir kinetic model to BTC data for ovalbumin (co = 1 g L-1) loaded on a Mustang XT5 module under binding conditions at flow rates of: 1.5 mL min-1 (■), 5 mL min-1 (▲), 10 mL min-1 (♦) and 20 mL min-1 (●).  the Mustang Q XT5 capsule such that a faster than realized approach to complete column loading is predicted (Figure 5.6). 5.4.4.2 Bi-Langmuir and SMA models As indicated by the sum of squared residuals (Table 5.3), use of either the bi-Langmuir (Equations 5-4 and 5-5) or the SMA (Equations 5-7 and 5-8) rate equations further improves the accuracy of the ZRM in matching experimental BTCs. For either of these more advanced rate equations, however, the improved fit is realized with a set of regressed model parameters that are physically unrealistic. To illustrate this point, the results for the bi-Langmuir implementation of the model will be discussed as it is the more accurate of the two. The ability to regress both the extra-column structural parameters of the ZRM and Dx to pulse or breakthrough data under flow-through conditions allows the parameters of the  103  Table 5.3  Sum of squared residuals for the fitting of various models to BTC data. Experimental conditions are the same as in Table 5.2. The performance of each binding model is that when implemented within the symmetric twozone ZRM.  Flow rate (mL min-1)  Langmuir  Bi-Langmuir  SMA  Spreading model  1.5  0.305  0.043  0.093  0.031  5  1.249  0.156  0.387  0.036  10  0.679  0.279  0.597  0.213  20  1.169  0.126  0.590  0.095  TOTAL  3.402  0.604  1.667  0.375  chosen binding rate model to be determined independently through global regression to a set of retained protein BTC data. The parameters for the bi-Langmuir model when implemented in the symmetric two-zone ZRM are shown in Table 5.4. Though the model represents the data set quite well (Figure 5.7), the two dynamic saturation binding capacities, q m ,1 and q m , 2 , change in unanticipated ways with flow rate: the capacity of higher affinity type 1 sites increases with flow rate while that of the lower affinity type 2 sites decreases, falling sharply between 10 and 20 mL min-1. A plausible explanation for this predicted extreme change in q m ,1 q m , 2 is not easily provided, as it suggests structural changes in membrane architecture that significantly alter the availabilities of the two different classes of binding sites. However, the membrane stack in the Mustang Q XT5 capsule is both mechanically supported and very thin (2.2 mm). As a result, the transmembrane pressure drop ∆P is extremely low, never exceeding 21 kPa (0.21 bar) over the range of flow rates recommended by the manufacturer and showing a linear dependence on interstitial velocity as predicted by the Blake-Kozeny equation [184]:  P  150  u int L 1   m   2  d p2 m  2  (5-23)  104  Table 5.4  Regressed parameters for the bi-Langmuir adsorption model when introduced into the two-zone ZRM. Experimental conditions and regression method are the same as in Table 5.2. Flow rate (mL min-1) 1.5 5 10 20  Figure 5.7  ka (L g-1 s-1)  kd (s-1)  ka,1 8.38x10-3  kd,1 3.18x10-5  ka,2 1.39x10-3  kd,2 5.49x10-4  q m ,1 (g L-1)  qm , 2 (g L-1)  152  195  181  184  224  171  304  43  Comparison of the symmetric two-zone ZRM incorporating the biLangmuir kinetic model to BTC data for ovalbumin (co = 1 g L-1) loaded on a Mustang XT5 module under binding conditions at flow rates of: 1.5 mL min-1 (■), 5 mL min-1 (▲), 10 mL min-1 (♦) and 20 mL min-1 (●).  105  which has been used previously to predict pressure drop across membranes [70] and represents a conservative estimate as it assumes spherical shaped voids with a surface area to volume ratio of av = 6/dp. Thus P is over-predicted by a factor of (3/2)2 for cylindrical pores where av = 4/dp. It is hard to reconcile the prediction of a large change in sorbent binding site composition with this result as such a change would presumably require gross changes to the architecture of the membrane that would be reflected in the dependence of ∆P on uint. For the model system under investigation here, this suggests that the use of the biLangmuir rate equation improves model performance not through an improved description of the intrinsic binding reaction, but rather through an increase in the number of regressed parameters (6 as opposed to 3). Similarly, globally regressed values for the capacity parameter of the SMA model ( qm       ) were found to increase from 317 mg mL-1 at 1.5 mL min-1 to 369 mg mL-1 at 20 mL min-1. As with the results for the biLangmuir isotherm, such a trend is not supported by the underlying model theory and suggests improved model accuracy is again simply due to a larger number of regressed parameters. 5.4.4.3 Spreading model When introduced into the symmetric two-zone ZRM, the spreading model provides a very accurate representation of the BTC data (Table 5.3). Model results match with experiment at all feed flow rates (Figure 5.8). Global regression of the general spreading model (Equations 5-9 and 5-10) to the data yields a simplified spreading model comprised of only two distinct adsorbed protein states.  Ovalbumin molecules are  predicted to adsorb reversibly to the sorbent (state 1) and then may reorient or spread reversibly to a second state (2) from which they cannot desorb directly into the bulk (i.e. setting ka,2 and kd,2 to zero results in no loss in model accuracy). The values of the regressed parameters (Table 5.5), which are consistent with those reported in a previous study [143] that applied a spreading-type model to IEXM chromatography, reveal that the simplified spreading model captures the relatively rapid increase in eluent concentration in the initial sections of each BTC by predicting relatively rapid population of state 1.  106  Figure 5.8  Comparison of the symmetric two-zone ZRM incorporating the simplified spreading model to BTC data for ovalbumin (co = 1 g L-1) loaded on a Mustang XT5 module under binding conditions at flow rates of: 1.5 mL min-1 (■), 5 mL min-1 (▲), 10 mL min-1 (♦) and 20 mL min-1 (●).  Table 5.5  Regressed parameters for the simplified spreading model when introduced into the two-zone ZRM. Experimental conditions and regression method are the same as in Table 5.2.  Flow rate (mL min-1)  ka,1  (L g-1 s-1)  kd,1  k12  (s-1)  (L g-1 s-1)  k21  K12  (s-1)  (M-1)  3.59x10-3  1.53x103  6.87x10-3  7.99x102  10  1.93x10-2  2.85x102  20  1.28x10-1  4.29x101  1.5 5  -3  8.28x10  -4  2.09x10  1.24x10  -4  qm ,1 (g L-1)    379 2.55  343  107  The pronounced asymmetry of BTCs at lower flow rates (1.5 and 5 mL min-1), a common observation in membrane chromatography systems, is related to the slower reaction kinetics associated with state 2, which can only become populated (and subsequently depopulated at high surface coverage, see below) following population of state 1. At higher feed flow rates (10 and 20 mL min-1), the BTCs become more symmetric, and the spreading model captures this through the flow-rate dependence of the reverse rate constant, k21, for the reorientation/spreading reaction.  In particular, while k12 is  insensitive to uint, k21 increases with uint. To understand this and connect it to observed breakthrough behavior, note that the spreading model computes the equilibrium amount of protein bound to a sorbent of known surface area by minimizing the Gibbs energy of the system for a given total protein concentration.  In the case of a non-  spherical/deformable protein that can bind in two distinct states, the amount of protein bound in each state at equilibrium therefore depends on the stability of each bound state, given by ∆G1 = - RT ln K1a for the lower affinity vertical/native state 1 and by ∆G2 = - RT ln K1aK12 for the higher affinity horizontal/spread state 2, normalized by the projected surface area per bound molecule (1 for state 1 and  = 2.55 for state 2). The difference between these two normalized energies is plotted in Figure 5.9 and shows that the weaker affinity state 1 will be more populated at breakthrough because it allows the proteins to load to a higher surface density.  Thus, as the surface approaches saturation, adsorbed  proteins in state 2 will revert to state 1 to free sorbent surface area and thereby permit the net energetically favorable loading of more protein in state 1. Figure 5.9 also shows that the overall driving force for conversion from state 2 back to state 1 increases with increasing uint. However, since an individual protein molecule bound in state 2 is at a lower energy than when bound in state 1, the rate of conversion back to state 1 depends on the availability of system energy that can be used by the molecule to overcome the affinity difference. One potential source of energy is that  provided by the wall shear stress  w , which for laminar flow in a pore generates a radially dependent force parallel to the pore wall.  By considering flow through a  membrane pore in this manner, one can show that:  108  Figure 5.9  Difference in the stability of bound protein states normalized by the projected surface area per bound molecule (G1 – G2/) for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 column at flow rates of 1.5, 5, 10 and 20 mL min-1.    w   where rp is the pore radius.  4 Q  mrp3  (5-24)   Thus,  w increases with Q, providing one plausible  explanation for the predicted increase in k21 with uint and the associated improvement in BTC symmetry. Indeed, this basic argument has been used to explain the influence of hydrodynamics on other biological binding equilibria at a fluid-solid interface [185]. These results suggest that a re-orientation of the adsorbed protein molecules is a much more plausible interpretation of the model than is a change in conformation, which would require a more substantial energy source than wall shear stress to reverse. More importantly, the spreading model implementation of the ZRM also provides an explanation for the common observation in IEXM chromatography of an increase in  109  dynamic binding capacity q* with increasing uint. Figure 5.10 reports experimental q* values for each of the BTC experiments performed and shows that q* for ovalbumin increases with increasing flow rate between 0 and 10 mL min-1. These BTC experiments were performed at a co = 1 g L-1, which is deep within the linear region of the static isotherm but does not correspond to saturation conditions.  Also shown in Figure 5.10  (dashed line) is the saturation capacity (338 g L-1) determined from (near-)static isotherm experiments, as well as model predictions of q*, q1, and q2. The model captures the curious increase in q* with flow rate, predicting that it occurs as a result of the depopulation of state 2 with increasing wall shear stress. It also yields a value for q m ,1 of 379 g L-1 which is constant and equal to qm,1 at Q = 10 mL min-1 and below. Thus, q m ,1 ,  Figure 5.10  Experimentally derived dynamic binding capacities, q*, (■) for the frontal loading of ovalbumin (co = 1 g L-1) on a Mustang XT5 column at flow rates of 1.5, 5, 10 and 20 mL min-1. Two-zone ZRM/spreading model predictions of q*, the concentrations of bound protein in state 1 (q1) and state 2 (q2), and the theoretical dynamic capacity based on all protein binding in state 1 ( q m ,1 ) are also shown. The dashed line represents the saturation capacity (qm) of the membrane for ovalbumin determined from (near-)static isotherm experiments.  110  which represents the hypothetical saturation capacity of the membrane if all ovalbumin bound in state 1, is higher than the measured static qm for the real system, which is predicted to include a small fraction of ovalbumin bound in state 2 (Figure 5.10). Increasing Q decreases q2, increasing both q1 and q*. Further increases in Q above 10 mL min-1 then lead to the expected decrease in dynamic binding capacity due to decreased retention time in the column. The spreading model therefore not only captures, but also provides a realistic underlying theory to explain, key experimental trends in the model system of ovalbumin loading onto a Mustang Q XT5 capsule. The spreading model predicts that the asymmetric nature of the observed BTCs is due to a series of processes at the sorbent surface, each characterized by unique kinetics. It has then been speculated that these processes may relate to changes in orientation or conformation of the protein upon adsorption, in accordance with recent experimental work that directly measures such changes with time for proteins adsorbed on an ionexchange resin [186,187]. However, the sequential surface kinetics defined by the model could relate to other physical phenomena such as surface aggregation. Finally, it is possible to now confirm the validity of omitting film mass transfer effects in the ZRM. Based on ka,1, Dkf is on the order of 5 x 104 across all flow rates, confirming that the rate of film mass transfer is orders of magnitude faster than the rate of solute adsorption.  5.2  Concluding remarks  It has been shown that the quantitative description and accurate prediction of IEXM chromatography can be achieved using the novel ZRM that properly accounts for regional contributions to band broadening so as to permit rigorous evaluation of putative binding models and associated binding mechanisms. For the loading of ovalbumin onto a Mustang Q XT5 IEXM capsule over a range of flow rates (1.5 - 20 mL min-1), this model-building strategy was used to characterize four potential solute binding models – Langmuir, bi-Langmuir, SMA and spreading. The spreading model, which describes a reversible change in orientation and/or conformation of protein molecules upon binding  111  to the pore surface, was found to accurately fit the asymmetric experimental breakthrough data and also provide mechanistic insight into the dependence of dynamic binding capacity on flow rate: a dependence that, although seen previously, has not been well explained. A spreading-type model similar to one proposed for ovalbumin binding to the XT5 IEXM capsule was recently used by Morbidelli and co-workers to describe experimentally observed elution patterns of human serum albumin from ion-exchange resin [188,189]. Thus, our ever increasing ability to more rapidly and accurately model complex transport processes occurring in preparative chromatography systems may be revealing the general utility of spreading-type binding models. However, though the spreading model was found most appropriate for describing the intrinsic binding process in this particular model membrane chromatography system, this will certainly not hold true for all systems or operational conditions. For example, the biLangmuir [131] or SMA [147] might indeed be the more appropriate binding model in certain cases. The value of this work is therefore in demonstrating the ability to use the ZRM to accurately assess the fundamental value of a putative binding model, as well as to isolate and define all contributions to band broadening occurring both within and external to the column. When used in the procedural framework outlined in this work, the ZRM therefore represents a flexible and useful theoretical tool to characterize membrane chromatography columns and to optimize and scale-up their operation so as to avoid unfavorable conditions while maintaining the intrinsic advantages of membrane chromatography.  112  Chapter 6  6 Conclusion As discussed in Chapter 1 considerable effort has been made to increase mammalian cell culture titres from levels in the low milligrams to their current concentrations of multigrams per liter [5]. These dramatic improvements will allow for > 100 kg batch sizes in the foreseeable future and have shifted production bottlenecks towards downstream processing where technological advances have not kept pace with upstream productivity increases [4,82].  There is therefore an increasing emphasis on the streamlining of  product recovery and purification in order to reduce costs in downstream processing which account for up to 80% of total manufacturing costs [15]. This streamlining can generally be accomplished by either reducing the total number of unit operations or by increasing their step yields. One method for reducing the number of unit operations is through rational process integration wherein two or more unit operations are combined into one step. CSAF is a novel technology that integrates the cell removal and product capture stages of the downstream processing of biopharmaceuticals by placing a specially designed rotor above a membrane chromatography column. The rotor assembly decouples TMP from the shear stress at the membrane surface consequently minimizing membrane fouling and achieving a uniform filtrate flux with no associated shear-induced cell lysis [16].  The  technology showed promising results at the bench-scale but was limited by a lack of fundamental phenomenological knowledge of its operation. The objective of this thesis therefore was to construct a framework of modeling tools that would allow for the refinement, scale-up and optimization of the CSAF design and operating parameters in silico. Firstly, a CFD model was built to describe the rotor chamber hydrodynamics of the laboratory-scale CSAF device. This model not only illuminated certain flaws in the original CSAF design – flaws that limited the degree to which the technology could be  113  scaled-up - but also provided a pathway along which a new optimal design could be developed: the rapid screening of multiple potential CSAF designs in silico. This CFD model was then used to scale-up CSAF, once again in silico, to an effective filtration area of ca. 6 m2. This represents an increase of more than 1,000 fold over what could have been achieved using the original design and permits the complete processing in ca. 2 h of 1000 L of culture harvested from either a perfusion, fed-batch or batch bioreactor. Once a modeling framework was constructed that would allow for the scale-up of the rotor assembly to an industrially relevant size it was important to develop a methodology to determine  the  optimum  operating  parameters  for  the  associated  membrane  chromatography column. To this end the novel ZRM was developed to describe solute mass transfer in stacked membrane chromatography.  This model expanded on the  current state of the art by explicitly accounting for all transport non-idealities within and external to the membrane stack including inhomogeneous radial flow distribution and system dispersion. Through its careful accounting of all major contributions to band broadening the ZRM was also shown to provide a useful framework for characterizing putative protein binding mechanisms and models for predicting BTCs and complex elution behavior, including the common observation that the dynamic binding capacity can increase with linear velocity in IEXM systems, and for simulating and scaling separations using membrane chromatography such as CSAF. By combining clarification with product capture and concentration in a single unit CSAF has the potential to greatly impact downstream processing of biopharmaceuticals. The technology could increase overall product yield by 10-20%, as well as increase the degree of automation and reduce processing volumes and operating and capital costs. The development in this thesis of a modeling framework to optimize and scale-up CSAF is an important and essential first step towards the implementation of the technology in an industrial manufacturing setting.  And in fact, due to the expense and complexity  involved in constructing multiple industrial-scale CSAF prototypes and then running multiple protein breakthrough experiments on them, modeling may represent the only feasible method of scale-up. However, modeling cannot fully replace data and there still exists a lack of experimental results essential in the development of CSAF.  114  For instance, although outside the scope and resources of this thesis project, experimentation with a full-scale CSAF unit is required to yield important information for the determination of an operational space for the technology. Of specific interest would be flux vs. TMP curves at different rotational velocities as well as studies investigating the effect of rotational frequency on cell viability. Furthermore, in order for the ZRM to effectively describe breakthrough from the large-scale CSAF unit, system dispersion and breakthrough curve data will have to be attained in order to be used in model building. Most importantly, an industrial-scale CSAF prototype will allow for a critical assessment of the technology consisting of a side-by-side comparison of the product yields and purities produced by the CSAF device to those produced by a sequence of conventional cell separation and product capture methods.  Such an  assessment will ensure that the CSAF represents a viable process alternative for the downstream processing of biopharmaceuticals. The ZRM also shows promise, as a stand-alone entity, to be an important tool for the characterization of any membrane chromatography unit that exhibits significant flow non-idealities.  Up to now the ZRM, in order to provide a sound basis for model  development and verification, has been tested on relatively simple systems: singlecomponent protein solutions, a flat sheet membrane stack, and conventional flow distribution systems. However, the ZRM could be generalized so as to be applicable over a wider range of systems, including industrial settings. Specifically such work should focus on expanding the model to multi-component mixtures and different flow geometries such as radial flow modules. This last area would be of particular interest since radial flow adsorbers are the most often encountered in industrial settings but current literature on membrane chromatography is focused almost entirely on flat sheet or hollow fibre membrane configurations [69]. Generalizing the ZRM to multi-component feedstocks processed through spiral wound membrane modules will therefore make the model more relevant to industry where such systems are the norm. The ZRM could also be expanded to describe traditional, bead-based chromatography. Although most models for separation in packed bed chromatography assume radially uniform axial velocities, flow distributions in chromatographic columns have been shown  115  to be non-uniform due to inadequate header design and/or non-uniform packing [168,170,190,191]. This radial heterogeneity in column properties and/or axial velocities would be well captured by the ZRM’s virtual radial zones. The ZRM could, therefore, represent a useful tool to describe band broadening in conventional chromatographic separations as well as membrane-based systems. The spreading model of protein adsorption developed in Chapter 5 also warrants further study. This model posited that the observed trends in the experimental BTCs were due to changes in conformation/orientation of the protein upon adsorption to the sorbent surface. This hypothesis could be experimentally tested using a variety of techniques. Massloading kinetics of protein adsorption can be measured using such highly sensitive methods as ellipsometry [192], surface plasmon resonance [193], total internal reflection fluorescence [192] and fluorescence recovery after photobleaching [194].  These  techniques are primarily used to determine adsorbed amounts but can also offer some clues concerning adsorption induced conformational changes [192]. Technologies also exist that have been previously used to more specifically characterize conformational changes of adsorbed proteins; these include circular dichroism spectroscopy [195], Raman spectroscopy [196], Fourier transform infrared spectroscopy [197] and fluorescence spectroscopy [192]. Differential scanning calorimetry can also provide insight into adsorbed protein energies and conformations [181] and has been used, as a complementary approach to the optical techniques mentioned above, to investigate adsorption-induced conformational changes of proteins [197]. Knowledge of surface energy changes would also improve understanding of the adsorption process. Although technically challenging, methods do exist for the experimental determination of surface energy during adsorption processes including axisymmetric drop shape analysis [198], measuring the radius of curvature of an atomic force microscopy micro-fabricated silicon nitride cantilever [199] and the recently developed resonant elastomeric surface-tension sensor technology [145]. Experimental observations of the adsorption process will allow for the refinement of the adsorption kinetic equations used as part of the ZRM and will therefore strengthen the model’s ability to quantitatively predict the performance of both stand-alone columns as well as a scaled-up CSAF unit.  116  One specific area where the implementation of CSAF technology shows particular promise is in the production of monoclonal antibodies (mAB). MABs are used in a variety of indications such as cancer, macular degeneration, multiple sclerosis, rheumatoid arthritis and Crohn’s disease [4,200] and as such represent a large and growing sector of biopharmaceutical production. In 2009 mAB-based products recorded global sales of $38 billion US and represented four of the top five selling biopharmaceuticals [4]. Furthermore the market for mABs has a large future potential: with approximately 240 antibodies or antibody derivatives currently in clinical trials such therapies dominate the biopharmaceutical pipeline [4,200,201] and their market is expected to rise to $56 billion US by 2012 [202]. Because antibodies are marketed for chronic conditions and have relatively low potency they typically require high dosages over a long period of time – 2-5 g per patient per year [11,203]. Consequently annual production levels can reach hundreds of kilograms per year, as opposed to the ~1 kg of EPO produced each year [7]. These large production levels make cost of manufacturing important especially considering that the cost of goods for mAB therapies can reach upwards of $1000 US per gram which makes treatment with these therapies prohibitively expensive.  Increased effort is therefore being made to  reduce manufacturing costs at the commercial scale by a factor of 10 or even 100 [203]. As previously discussed, one way to reduce manufacturing costs and thus meet market demand at an acceptable cost to the consumer is through the type of unit integration seen in CSAF technology. The membrane chromatography column within CSAF makes it especially well-suited for mAB production. Although traditional bead based chromatography is the predominant paradigm in mAB purification much work is being done to investigate membrane-based alternatives (see Boi’s review article [204]). Above and beyond the mass transfer issues discussed throughout this thesis membranes offer a larger available surface area than traditional chromatographic beads for binding large protein molecules such as mABs because these large molecules cannot enter the small pores of the beads [204]. Furthermore, the membranes within CSAF give it a desirable flexibility. The membranes can either be of an affinity modality (e.g. Protein A) - currently the most widely used  117  methodology for mAB purification [205-207] - or of an ion-exchange modality which is increasingly being investigated as a capture step due to the disadvantages associated with Protein A (e.g. leaking and high costs) [83]. These advantages give CSAF technology the potential for being an integral tool in mAB production.  118  References [1]  Cohen S, Chang CY, Boyer HW, Helling RB. 1973. Construction of biologically functional bacterial plasmids in vitro. Proc Natl Acad Sci 70:3240-3244.  [2]  Kohler G, Milstein C. 1977. Continuous cultures of fused cells secreting antibody of predefined specificity. Nature 256:495-497.  [3]  Walsh, G. 2003. Pharmaceutical biotechnology products approved within the European Union. Eur J Pharm Biopharm 55:3-10.  [4]  Walsh G. 2010. Biopharmaceutical benchmarks 2010. Nat Biotechnol 28:917924.  [5]  Wurm F. 2004. Production of recombinant protein therapeutics in cultivated mammalian cells. Nat Biotechnol 22:1393-1398.  [6]  Chu L, Robinson DK. 2001. Industrial choices for protein production by largescale cell culture. Curr Opin Biotechnol 12:180-187.  [7]  Kretzmer G. 2002. Industrial processes with animal cells. Appl Microbiol Biotechnol 59:135-142  [8]  Warnock JN, Al-Rubeai M. 2006. Bioreactor systems for the production of biopharmaceuticals from animal cells. Biotechnol Appl Biochem 45:1-12.  [9]  Castilho LR, Medronho RA. 2002. Cell retention for suspended-cell perfusion cultures. Adv Biochem Eng 74:129-169.  [10]  Butler M. 2005. Animal cell cultures: recent achievements and perspectives in the production of biopharmaceuticals. App Microbiol Biotechnol 65: 283-291.  [11]  Sommerfeld S, Strube J. 2005. Challenges in biotechnology production – generic processes and process optimization for monoclonal antibodies. Chem Eng Process 44: 1123-1137.  [12]  Lubinecki AS. 1998. Historical reflections on cell culture engineering. Cytotechnology 28:139-145.  [13]  Jagchies G, Gronberg A, Bjorkman T, Lacki K. Johansson HJ. 2006. Technical and economical evaluation of downstream processing options for monoclonal antibody (Mab) production. Biopharm Int10-15.  [14]  Low D. 2005. Process Development: Time to Innovate or Time to Market? Crossroads of Biotechnology Conference. Montreal, QC. 119  [15]  Lowe CR, Lowe G, Gupta G. 2001. New developments in affinity chromatography with potential application in the production of biopharmaceuticals. J Biochem Biophys Methods 49:561-574.  [16]  Vogel JH, Anspach B, Kroner K, Piret JM, Haynes CA. 2002. Controlled shear affinity filtration (CSAF): a new technology for integration of cell separation and protein isolation from mammalian cell cultures. Biotechnol Bioeng 78:806-814.  [17]  Batt B, Davis RH, Kompala DS. 1990. Inclined sedimentation for selective retention of viable hybridomas in a continuous suspension bioreactor. Biotechnol Prog 6:458-464.  [18]  Voisard D, Meuwly F, Ruffieux P-A, Baer G, Kadouri A. 2003. Potential of cell retention techniques for large-scale high-density perfusion culture of suspended mammalian cells. Biotechnol Bioeng 82:751-765.  [19]  Woodside SM, Bowen BD, Piret JM. 1998. Mammalian cell retention devices for stirred perfusion bioreactors. Cytotechnology 28:163-175.  [20]  Trampler F, Sonderhoff SA, Pui PWS, Kilburn DG, Piret JM. 1994. Acoustic cell filter for high density perfusion culture of hybridoma cells. Biotechnol 12:281284.  [21]  Shirgaonkar IZ, Lanthier S, Kamen A. 2004. Acoustic cell filter: a proven cell retention technology for perfusion of animal cell cultures. Biotechnol Adv 22:433-444.  [22]  Deo YM, Mahadevan MD, Fuchs R. 1996. Practical considerations in operation and scale-up of spin filter based bioreactors for monoclonal antibody production. Biotechnol Prog 12:57-64.  [23]  Himmelfarb P, Thayer PS, Martin HE. 1969. Spin-filter culture: the propagation of mammalian cells in suspension. Science 164:555-557.  [24]  Iding K, Lütkemeyer D, Fraune E, Gerlach K, Lehmann J. 2000. Influence of alterations in culture condition and changes in perfusion parameters on the retention performance of a 20 m spinfilter during a perfusion cultivation of a recombinant CHO cell line in pilot scale. Cytotechnology 34:141-150.  [25]  Aunins JG, Wang IC. 1989. Induced flocculation of animal cells in suspension. Biotechnol Bioeng 34:629-638.  [26]  de la Broise D, Noiseux M, Lemieux R. 1991. Long-term perfusion culture of hybridoma: a “grow or die” cell cycle system. Biotechnol Bioeng 38:781-787.  120  [27]  Maiorella B, Dorin G, Carion A, Harano D. 1991. Crossflow microfiltration of animal cells. Biotechnol Bioeng 37:121-126.  [28]  van Reis R, Leonard LC, Hsu CC, Builder SE. 1991. Industrial scale harvest of proteins from mammalian cell culture by tangential flow filtration. Biotechnol Bioeng 38:413-422.  [29]  Zhang S, Handa-Corrigan A, Spier RE. 1993. A comparison of oxygenation methods for high-density perfusion cultures of animal cells. Biotechnol Bioeng 41:685-692.  [30]  Takazawa Y, Tokashiki M. 1898. High cell density perfusion culture of mousehuman hybridoma. Appl Microbiol Biotechnol 32:280-284.  [31]  Maia ABRA, Nelson DL. 1993. Application of gravitational sedimentation to efficient cellular recycling in continuous culture with cell recycle. Biotechnol Bioeng 41:361-369.  [32]  Su WW, He BJ, Liang H, Sun S. 1996. A perfusion airlift bioreactor for highdensity plant-cell cultivation and secreted protein-production. J Biotechnol 50:225-233.  [33]  Hülscher M, Scheibler U, Onken U. 1992. Selective recycle of viable animal cells by coupling of airlift reactor and cell settler. Biotechnol Bioeng 39:442-446.  [34]  Arai T, Yokohama S, Tokashiki M. 1993. 50 l scale perfusion culture of hybridoma cells by gravitational settling for cell separation. In: Kaminogawa S, Ametani A, Hachimura S, editors. Animal cell technology: basic and applied aspects. Dordnecht: Kluwer Academic Publishers. p. 341-346.  [35]  Wen ZY, Teng XW, Chen F. 2000. A novel perfusion system for animal cell cultures by two step sequential sedimentation. J Biotechnol 79:1-11.  [36]  Searles JA, Todd P, Kompala DS. 1994. Viable cell recycle with an inclined settler in the perfusion culture of suspended recombinant Chinese hamster ovary cells. Biotechnol Prog 10:198-206.  [37]  Knaack C, André G, Chavarie C. 1994. Conical bioreactor with internal lamella settler for perfusion culture of suspension cells. In: Spier RE, editor. Animal cell technology. Oxford: Butterworth-Heinemann. p. 230-233.  [38]  Thompson KJ, Wilson JS. 1994. A compact gravitational settling device for cell retention. In: Spier RE, editor. Animal cell technology. Oxford: ButterworthHeinemann. p. 227-229.  121  [39]  Stevens J, Eickel S, Onken U. 1994. Lamellar clarifier – a new device for animal cell retention in perfusion culture. In: Spier RE, editor. Animal cell technology. Oxford: Butterworth-Heinemann. p. 234-239.  [40]  Erikson RA. 1984. Disk stack centrifuges in biotechnology. Chem Eng Prog 80:51-54.  [41]  Tokashiki M, Arai T, Hamamoto K, Ishimaru K. 1990. High density culture of hybridoma cells using a perfusion culture vessel with an external centrifuge. Cytotechnology 3:239-244.  [42]  Jäger V. 1992. A novel perfusion system for the large-scale cultivation of animal cells based on a continuous flow centrifuge. In: Spier RE, Griffiths JB, MacDonald C, editors. Animal cell technology: Developments, processes, products. Oxford: Butterworth-Heinemann. p. 397-402.  [43]  Al-Rubeai M, Singh RP, Goldman MH, Emery AN. 1995. Death mechanisms of animal cells in conditions of intensive agitation. Biotechnol Bioeng 45:465-472.  [44]  Keane J, Ryan D, Gray P. 2003. Effect of shear stress on expression of a recombinant protein by Chinese hamster ovary cells. Biotechnol Bioeng 81:211220.  [45]  Kretzmer G, Schurgel K. 1991. Response of mammalian cells to shear stress. Appl Microbiol Biotechnol 34:613-616.  [46]  Ludwig A, Kretzmer G. 1993. Shear stress induced variation of cell condition and productivity. J Biotechnol 27:217-223.  [47]  Mardikar S, Niranjan K. 2000. Observations on the shear damage to different animal cells in a concentric cylinder viscometer. Biotechnol Bioeng 68:697-704.  [48]  McQueen A, Meilhoc E, Bailey J. 1987. Flow effects on the viability and lysis of suspended mammalian cells. Biotechnol Lett 9:831-836.  [49]  van der Pol L, Tramper J. 1998. Shear sensitivity of animal cells from a culturemedium perspective. Trends Biotechnol 16:323-328.  [50]  Yavorsky D, Blanck R, Lambalot C, Brunkow R. 2003. The clarification of bioreactor cell cultures for biopharmaceuticals. Pharma Technol 27:62-76.  [51]  van Reis R, Zydney A. 2001. Membrane separations in biotechnology. Curr Op Biotechnol 12:208-211.  [52]  Russell E, Wang A, Rathore AS. 2007. Harvest of a therapeutic protein product from high cell density fermentation broths: principles and case study. In: Shukla  122  AA, Etzel MR, Gadam S, editors. Process scale bioseparations for the biopharmaceutical industry. Boca Raton: CRC Press. p. 1-58. [53]  Kawahara H, Mitsuda S, Kumazawa E, Takeshita Y. 1994. High density culture of FM-3a cells using a bioreactor with an external tangential flow filtration device. Cytotechnology 24:61-66.  [54]  Vogel JH, Kroner K. 1999. Controlled shear filtration: a novel technique for animal cell separation. Biotechnol Bioeng 63:663-674.  [55]  Zijlstra GM, Michielsen MJF, de Gooijer CD, van der Pol LA, Tramper J. 1998. IgG and hybridoma partitioning in aqueous two-phase systems containing a dyeligand. Bioseparation 7:117-126.  [56]  Selber K, Tjerneld F, Collen A, Hyytia T, Nakari-Setala T, Bailey M, Fagerstrom R, Kan J, van der Laan J, Penttila M, Kula M-R. 2004 Large-scale separation and production of engineered proteins, designed for facilitated recovery in detergentbased aqueous two-phase extraction systems. Proc Biochem 39:889-896.  [57]  Rito-Palomares M. 2004. Practical application of aqueous two-phase partition to process development for the recovery of biological products. J Chromatogr B 807:3-11.  [58]  Benavides J, Aguilar O, Lapizco-Encinas BH, Rito-Palomares M. 2008. Extraction and purification of bioproducts and nanoparticles using aqueous twophase systems strategies. Chem Eng Technol 31:838-845.  [89]  Rosa PAJ, Azevedo AM, Sommerfeld S, Mutter M, Aires-Barros MR, Bäcker W. 2009. Application of aqueous two-phase systems to antibody purification: a multistage approach. J Biotechnol 139:306-313.  [60]  Dainiak MB, Izumrudov VA, Muronetz VI, Galaev IY, Mattiason B. 1999. Affinity precipitation of monoclonal antibodies by nonstoichiometric polyelectrolyte complexes. Bioseparation 7:231-240.  [61]  Taipa MA, Kaul R, Mattiason B, Cabral JMS. 1998. Preliminary studies on the purification of a monoclonal antibody by affinity precipitation with Eudragit S100. J Molec Recog 11:240-242.  [62]  Lim S, Manusu HP, Gooley AA, Williams KL, Rylatt DB. 1998. Purification of monoclonal antibodies from ascitic fluid using preparative electrophoresis. J Chromatogr A 827:329-335.  [63]  Thomas TM, Shave EE, Bate IM, Gee SC, Franklin S, Rylatt DB. 2002. Preparative electrophoresis: a general method for the purification of polyclonal antibodies. J Chromatogr A 944:161-168.  123  [64]  Lee HG. 1997. Rapid high-performance isoelectric focusing of monoclonal antibodies in uncoated fused-silica capillaries. J Chromatogr A 790:215-223.  [65]  Naveh D, Siegel RC. 1991. Large scale downstream processing of monoclonal antibodies. Bioseparation 1:351-366.  [66]  Thömmes J, Kula M-R. 1995. Membrane chromatography – an integrative concept in the downstream processing of proteins. Biotechnol Prog 11:357-367.  [67]  van Reis R, Zydney A. 2007. Bioprocess membrane technology. J Mem Sci 297:16-50.  [68]  Klein E. 2000. Affinity membranes: a 10-year review. J Mem Sci 179:1-27.  [69]  Ghosh R. 2002. Protein separation using membrane opportunities and challenges. J Chromatogr A 952:13-27.  [70]  Suen S-Y, Etzel M. 1992. A mathematical analysis of affinity membrane bioseparations. Chem Eng Sci 47:1355-1364.  [71]  Charcosset C. 2006 Membrane processes in biotechnology: an overview. Biotechnol Adv 24:482-492.  [72]  Niazi SK. 2006. Handbook of biogeneric therapeutic proteins. Boca Raton: CRC Press. p. 217-236.  [73]  Vailaya A. 2005. Fundamentals of reversed-phase chromatography: thermodynamic and extrathermodynamic treatment. J Liq Chromatogr Relat Technol 28:965-1054.  [74]  Shukla AA, Yigzaw Y. 2007. Modes of preparative chromatography. In: Shukla AA, Etzel MR, Gadam S, editors. Process scale bioseparations for the biopharmaceutical industry. Boca Raton: CRC Press. p. 179-226.  [75]  Brandt S, Goffe RA, Kessler SB, O’Connor JL, Zale SE. 1988. Membrane based affinity technology for commercial scale purifications. Biotechnol (NY) 6:779782.  [76]  Saxena A, Tripathi BP, Kumar M, Shahi VK. 2009. Membrane-based techniques for the separation and purification of proteins: an overview. Adv Colloid Interface Sci 145:1-22.  [77]  Castilho LR, Anspach FB, Deckwer W. 2002. Comparison of affinity membranes for the purification of immunoglobulins. J Mem Sci 207:253-264.  chromatography:  124  [78]  Roper DK, Lightfoot EN. 1995. Separation of biomolecules using adsorptive membranes. J Chromatogr A 702:3-26.  [79]  Zhou JX, Tressel T, Gottschalk U, Solamo F, Pastor A, Dermawan S, Hong T, Reif O, Mora J, Hutchinson F, Murphy M. 2006. New Q membrane scale-down model for process-scale antibody purification. J Chromatogr A 1134:66-73.  [80]  Liu H-C, Fried JR. 1994. Breakthrough of lysozyme through an affinity membrane of cellulose-cibracon blue. AIChE J 40:40-49.  [81]  Zeng X, Ruckenstein E. 1999. Membrane chromatography: preparation and applications to protein separation. Biotechnol Prog 15:1003-1019.  [82]  Gottschalk U. 2008. Bioseparation in antibody manufacturing: the good, the bad and the ugly. Biotechnol Prog 24:496-503.  [83]  Low D, O’Leary R, Pujar NS. 2007. Future of antibody purification. J Chromatogr B 848:48-63.  [84]  Hubbuch J, Thommes J, Kula M-R. 2005. Biochemical engineering aspects of expanded bed adsorption. Adv Biochem Eng/Biotechnol 92:101-123.  [85]  Hart RA, Lester PM, Reifsnyder DH, Ogez JR, Builder SE. 1994. Large scale, in situ isolation of periplasmic IGF-I from E. coli. Biotechnol 12:1113-1117.  [86]  Persson J, Andersen DC, Lester PM. 2005. Evaluation of different primary recovery methods for E. coli-derived recombinant human growth hormone and compatibility with further down-stream processing. Biotechnol Bioeng 90:442451.  [87]  Przybycien TM, Pujar NS, Steele LM. 2004. Alternative bioseparation operations: life beyond packed-bed chromatography. Curr Op Biotechnol 15:469-478.  [88]  Anspach FB, Curbelo D, Hartmann R, Garke G, Deckwer W-D. 1999. Expandedbed chromatography in primary protein purification. J Chromatogr A 865:129144.  [89]  Chase H. 1994. Purification of proteins by adsorption chromatography in expanded beds. Trends Biotechnol 12:296-303.  [90]  Fahrner RL, Blank GS, Zapata GA. 1999. Expanded bed protein A affinity chromatography of a recombinant humanized monoclonal antibody: process development, operation and comparison with a packed bed method. J Biotechnol 75:273-280.  125  [91]  Ohashi R, Otero JM, Chwistek A, Yamato I, Hamel J-PP. 2002. On-line purification of monoclonal antibodies using an integrated stirred-tank reactor/expanded bed adsorption system. Biotechnol Prog 18:1292-1300.  [92]  Thömmes J, Bader A, Halfar M, Karau A, Kula M-R. 1996. Isolation of monoclonal antibodies from cell containing hybridoma broth using a protein A coated adsorbent in expanded beds. J Chromatogr A 752:111-122.  [93]  Niazi SK. 2006. Handbook of biogeneric therapeutic proteins. Boca Raton: CRC Press. p. 189-216.  [94]  Valdes R, Gomez L, Padilla S, Brito J, Reyes B, Alvarez T, Mendoza O, Herrera O, Ferro W, Pujol M, Leal V, Linares M, Hevia Y, Garcia C, Mila L, Garcia O, Sanchez R, Acosta A, Geada D, Paez R, Vega JL, Borroto C. 2003. Large-scale purification of an antibody directed against hepatitis B surface antigen from transgenic tobacco plant. Biochem Biophys Res Commun 308: 94-100.  [95]  Lullau E, Kanttinen A, Hassel J, Berg M, Haag-Alvarsson A, Cederbrant K, Greenberg B, Fenge C, Schweikart F. 2003 Comparison of batch and perfusion culture in combination with pilot scale expanded bed purification for production of soluble recombinant -secretase. Biotechnol Prog 19:37-44.  [96]  Blank GS, Zapata G, Fahrner R, Milton M, Yedinak C, Knudsen H, Schmelzer C. 2001. Expanded bed adsorption in the purification of monoclonal antibodies: a comparison of process alternatives. Bioseparation 10:65-71.  [97]  Hubbuch J, Thomas ORT. 2002 High-gradient magnetic affinity separation of trypsin from porcine pancreatin. Biotechnol Bioeng 79:302-312.  [98]  Heeboll-Nielsen A, Justesen SFL, Hobley TJ, Thomas ORT. 2004. Superparamagnetic adsorbents for high-gradient magnetic fishing of lectins out of legume extracts. Biotechnol Bioeng 87:311-323.  [99]  Heeboll-Nielsen A, Justesen SFL, Thomas ORT. 2004. Fractionation of whey proteins with high-capacity superparamagnetic ion-exchanges. J Biotechnol 113:247-262.  [100] Unger DR, Muzzio FJ, Aunins JG, Singhvi R. 2000. Computational and experimental investigation of flow and fluid mixing in the roller bottle bioreactor. Biotechnol Bioeng 70:117-130. [101] Williams KA, Saini S, Wick TM. 2002. Computational fluid dynamics modeling of steady-state momentum and mass transport in a bioreactor for cartilage tissue engineering. Biotechnol Prog 18:951-963.  126  [102] Davidson KM, Sushil S, Eggleton CD, Marten MR. 2003. Using computational fluid dynamics software to estimate circulation time distributions in bioreactors. Biotechnol Prog 19:1480-1486. [103] Haut B, Ben Amor H, Coulon L, Jacquet A, Halloin V. 2003. Hydrodynamics and mass transfer in a Couette-Taylor bioreactor for the culture of animal cells. Chem Eng Sci 58:777-784. [104] Brannock MWD, De Wever H, Wang Y, Leslie G. 2009. Computational fluid dynamics simulations of MBRs: inside submerged versus outside submerged membranes. Desalination 236:244-251. [105] Mousavi SM, Jafari A, Yaghmaei S, Vossoughi M, Turunen I. 2008. Experiments and CFD simulation of ferrous biooxidation in a bubble column bioreactor. Comp Chem Eng 32:1681-1688. [106] Vial C, Poncin S, Wild G, Midoux N. 2002. Experimental and theoretical analysis of the hydrodynamics in the riser of an external loop airlift reactor. Chem Eng Sci 57:4745-4762. [107] Boychyn M, Yim SSS, Ayazi Shamlou P, Bulmer M, More J, Hoare M. 2001. Characterization of flow intensity in continuous centrifuges for the development of laboratory mimics. Chem Eng Sci 56:4759-4770. [108] Boychyn M, Yim SSS, Bulmer M, More J, Bracewell DG, Hoare M. 2004. Performance prediction of industrial centrifuges using scale-down models. Bioprocess Biosyst Eng 26:385-391. [109] Taha T, Cui ZF. 2002. CFD modeling of gas-sparged ultrafiltration in tubular membranes. J Mem Sci 210:13-27. [110] Tarabara VV, Wiesner MR. 2003. Computational fluid dynamics modeling of the flow in a laboratory membrane filtration cell operated at low recoveries. Chem Eng Sci 58:239-246. [111] Wu YX, Yu HW, Ching CB. 2004. A computational fluid dynamics study of binary adsorption separation in chromatography. Chem Eng Technol 27:955-961. [112] Bouzerar R, Jaffrin MY, Ding L, Paullier P. 2000. Influence of geometry and angular velocity on performance of a rotating disk filter. AIChE J 46:257-265. [113] Rainer M, Hoflinger W, Koch W, Pongratz E, Oechsle D. 2002. 3D-flow simulation and optimization of a cross flow filtration with rotating discs. Sep Pur Technol 26:121-131.  127  [114] Serra CA, Wiesner MR, Laine J-M. 1999. Rotating membrane disk filters: design evaluation using computational fluid dynamics. Chem Eng J 72:1-17. [115] Serra CA, Wiesner MR. 2000. A comparison of rotating and stationary membrane disk filters using computational fluid dynamics. J Mem Sci 165:19-29. [116] Castilho LR, Anspach FB. 2003. CFD-aided design of a dynamic filter for mammalian cell separation. Biotechnol Bioeng 83:515-524. [117] Williams M, Chen WC, Bache G, Eastland A. 1991. An analysis methodology for internal swirling flow systems with a rotating wall. J Turbomachinery 113:83-90. [118] Fluent Inc. 1998. Computational fluid dynamics software user’s guide. [119] Pope SB. 2000. Turbulent Flows. Cambridge: Cambridge University Press. p. 93. [120] Brinkman HC. 1942. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Appl Sci Res A 1: 27-34. [121] Hill DF, Sharp KV, Adrian RJ. 2000. Stereoscopic particle image velocimetry measurements of the flow around a Rushton turbine. Exp Fluids 29:478-485. [122] Hopkins LM, Kelly JT, Wexler AS, Prasad AK. 2000. Particle image velocimetry measurements in complex geometries. Exp Fluids 29:91-95. [123] Pruvost J, Legrand J, Legentilhomme P, Doubliez L. 2000. Particle image velocimetry investigation of the flow-field of a 3D turbulent annular swirling decaying flow induced by means of a tangential inlet. Exp Fluids 29:291-301. [124] Shafiqul Islam M, Haga K, Kaminaga M, Hine R, Monde M. 2002. Experimental analysis of turbulent flow structure in a fully developed rib-roughened rectangular channel with PIV. Exp Fluids 33:296-306. [125] Xiong W, Kalkuhler K, Merzkirch W. 2003. Velocity and turbulence measurements downstream of flow conditioners. Flow Meas Inst 14:249-260. [126] Armenante PM, Luo C, Chou C-C, Fort I, Medek J. 1997. Velocity profiles in a closed, unbaffled vessel: comparison between experimental LDV data and numerical CFD predictions. Chem Eng Sci 52:3483-3492. [127] Khopkar AR, Aubin J, Xuereb C, Le Sauze N, Bertrand J, Ranade VV. 2003. Gas-liquid flow generated by a pitched-blade turbine: particle image velocimetry measurements and computational fluid dynamics simulations. Ind Eng Chem Res 42:5318-5332.  128  [128] Ranade VV. 1997. An efficient computational model for simulating flow in stirred vessels: a case of Rushton turbine. Chem Eng Sci 52:4473-4484. [129] Sheng J, Meng H, Fox RO. 1998. Validation of CFD simulations of a stirred tank using particle image velocimetry data. Can J Chem Eng 76:611-625. [130] von Karman T. 1921. On laminar and turbulent friction. Zeits F angew Math u Mech 1:233-252. [131] Boi C, Dimartino S, Sarti GC. 2007. Modelling and simulation of affinity membrane adsorption. J Chromatogr A 1162:24-33. [132] Ghosh R, Wong T. 2006. Effect of module design on the efficiency of membrane chromatographic separation processes. J Mem Sci 281:532-540. [133] Sarfert FT, Etzel MR. 1997. Mass transfer limitations in protein separations using ion-exchange membranes. J Chromatogr A 764:3-20. [134] Briefs K-G, Kula M-R. 1992. Fast protein chromatography on analytical and preparative scale using modified microporous membranes. Chem Eng Sci 47:141149. [135] Gebauer KH, Thömmes J, Kula M-R. 1997. Breakthrough performance of highcapacity membrane adsorbers in protein chromatography. Chem Eng Sci 52:405419. [136] Roper DK, Lightfoot EN. 1995. Estimating plate heights in stacked-membrane chromatography by flow reversal. J Chromatogr A 702:69–80. [137] Shi W, Zhang F, Zhang G. 2005. Mathematical analysis of affinity membrane chromatography. J Chromatogr A 1081:156-162. [138] Suen S-Y, Etzel M. 1994. Sorption kinetics and breakthrough curves for pepsin and chymosin using pepstatin A affinity membranes. J Chromatogr A 686:179192. [139] Teeters MA, Root TW, Lightfoot EN. 2002. Performance and scale-up of adsorptive membrane chromatography. J Chromatogr A 944:129-139. [140] Vicente T, Sousa MFQ, Peixoto C, Mota JPB, Alves PM, Carrondo MJT. 2008. Anion-exchange membrane chromatography for purification of rotavirus-like particles. J Mem Sci 311:270–283. [141] Wang J, Dismer F, Hubbuch J, Ulbricht M. 2008. Detailed analysis of membrane adsorber pore structure and protein binding by advanced microscopy. J Mem Sci 320:456–467.  129  [142] Yang H, Bitzer M, Etzel MR. 1999. Analysis of protein purification using ionexchange membranes. Ind Eng Chem Res 38:4044-4050. [143] Yang H, Etzel MR. 2003. Evaluation of three kinetic equations in models of protein purification using ion-exchange membranes. Ind Eng Chem Res 42:890– 896. [144] Talbot J, Tarjus G, van Tassel PR, Viot P. 2000. From car parking to protein adsorption: an overview of sequential adsorption processes. Colloids Surfaces 165: 287-324. [145] Clark AJ, Kotlicki A, Haynes CA, Whitehead LA. 2007. A new model of protein adsorption kinetics derived from simultaneous measurement of mass loading and changes in surface energy. Langmuir 23:5591-5600. [146] Lundstrom I. Models of protein adsorption on solid surfaces. Prog Colloid Polymer Sci 70:76-82. [147] Brooks CA, Cramer SM. 1992. Steric mass-action ion exchange: displacement profiles and induced salt gradients. AIChE J 38:1969-1978. [148] Lightfoot EN, Athalye AM, Coffman JL, Roper DK, Root TW. 1995. Nuclear magnetic resonance and the design of chromatographic separations. J Chromatogr A 707:45–55. [149] Montesinos-Cisneros RM, de la Vega Olivas J, Ortega J, Guzmán R, TejedaMansir A. 2007. Breakthrough performance of plasmid DNA on ion-exchange membrane columns. Biotechnol Prog 23:881–887. [150] Montesinos-Cisneros RM, Lucero-Acuna A, Ortega J, Guzman R, Tejeda-Mansir A. 2007. Breakthrough performance of large proteins on ion-exchange membrane columns. Biotechnol Appl Biochem 48:117-125. [151] Dancette OP, Taboureau J-L, Tournier E, Charcosset C, Blond P. 1999. Purification of immunoglobulin G by Protein A/G affinity membrane chromatograph. J Chromatogr A 723:61-68. [152] Serafica GC, Pimbley J, Belfort G. 1994. Protein fractionation using fast flow immobilized metal chelate affinity membranes. Biotechnol Bioeng 43:21-36. [153] Shiosaki A, Goto M, Hirose T. 1994. Frontal analysis of protein adsorption on a membrane adsorber. J Chromatogr A 679:1-9  130  [154] Kim M, Saito K, Furusaki S. 1991. Adsorption and elution of bovine g-globulin using affinity membrane containing hydrophobic amino acids as ligands. J Chromatogr 585:45-51. [155] Saffman PG. 1964. The lift force on a small sphere in a slow shear flow. J Fluid Mech 22:385-400. [156] Belfort G, Davis RH, Zydney AL. 1994. The behavior of suspensions and macromolecular solutions in crossflow microfiltration. J Mem Sci 96:1-58. [157] Savins JG, Metzner AB. 1970. Radial (secondary) flows in a rheogoniometric devices. Rheol Acta 9:365-373. [158] Ellenberger J, Fortuin JMH. 1985. A criterion for purely tangential laminar flow in the cone-and-plate rheometer and the parallel-plate rheometer. Chem Eng Sci 40:111-116. [159] Patankar SV. 1980. Numerical heat transfer and fluid flow. Washington: Hemisphere Publishing Co. [160] Holland CD, Liapis AI. 1983. Computer Methods for Solving Dynamic Separation Problems. New York: McGraw-Hill Book Co. p. 17-19. [161] Sdougos HP, Bussolari SR, Dewey CF. 1984. Secondary flow and turbulence in a cone-and-plate device. J Fluid Mech 138:379-404. [162] Gill PE, Murray W. 1978. Algorithms for the solution of the nonlinear leastsquares problem. SIAM J Numer Anal 15:977-992. [163] Daily JW, Nece RE. 1960. Chamber dimension effects on induced flow and frictional resistance of enclosed rotating disks. ASME J Basic Engineering 82:217-232. [164] Ketola HN, McGrew JM. 1968. Pressure, frictional resistance, and flow characteristics of the partially wetted rotating disk. Trans ASME, J Lubrication Technol 90:395-404. [165] Avramescu M-E, Borneman Z, Wessling M. 2003. Dynamic behavior of adsorber membranes for protein recovery. Biotechnol Bioeng 84;564-572. [166] von Lieres E, Wang J, Ulbricht M. 2010. Model based quantification of internal flow distributions from breakthrough curves of flat membrane chromatography modules. Chem Eng Technol 33:960-968. [167] Barber JA, Perkins JD, Sargent RWH. 1998. Boundary conditions for flow with dispersion. Chem Eng Sci 53:1463–1464.  131  [168] Farkas T, Guiochon G. 1997. Contribution of the radial distribution of the flow veolocity to band broadening in HPLC columns. Anal Chem 69:4592-4600. [169] Shalliker RA, Broyles BS, Guiochon G. 2000. Physical evidence of two wall effects in liquid chromatography. J Chromatogr A 888:1-12. [170] Yun T, Guiochon G. 1997. Visualization of the heterogeneity of column beds. J Chromatogr A 760:17-24. [171] Schiesser WE. 1991. The Numerical Method of Lines. San Diego: Academic Press. [172] Shampine LF, Reichelt MW. 1997. The MATLAB ODE Suite. J SIAM Sci Comput 18:1-22. [173] Coleman TF, Li Y. 1996. An interior, trust region approach for nonlinear minimization subject to bounds. SIAM J Optim 6:418–445. [174] Coleman TF, Li Y. 1994. On the convergence of reflective Newton methods for large-scale nonlinear minimization subject to bounds. Math Program 67:189-224. [175] Tyn MT, Gusek TW. 1990. Prediction of diffusion coefficients of proteins. Biotechnol Bioeng 35:327-338. [176] De Ligny CL. 1970. Coupling between diffusion and convection in radial dispersion of matter by fluid flow through packed beds. Chem Eng Sci 47:253262. [177] Kochan JE, Wu Y-J, Etzel MR. 1996. Purification of bovine immunoglobulin G via protein G affinity membranes. Ind Eng Chem Res 35:150-1155. [178] Gritti F, Piatkowski W, Guiochon G. 2003. Study of the mass transfer kinetics in a monolithic column. J Chromatogr A 983:51-71. [179] Tejeda-Mansir A, Montesinos RM, Guzmán R. 2001. Mathematical analysis of frontal affinity chromatography in particle and membrane configurations. J Biochem Biophys Methods 49:1–28. [180] Daly SM, Przybycien TM, Tilton RD. 2003. Coverage-dependent orientation of lysozyme adsorbed on silica. Langmuir 19:3848-3857. [181] Haynes CA, Norde W. 1995. Structures and stabilities of adsorbed proteins. J Colloid Interface Sci 169:313-328.  132  [182] Sethuraman A, Vedantham G, Imoto T, Przybycien TM, Belfort G. 2004. Protein unfolding at interfaces: Slow dynamics of -helixto b-sheet transition. Proteins: Struct, Funct, Bioinf 56:669-678. [183] Athalye AM, Gibbs SJ, Lightfoot EN. 1991. Predictability of chromatographic protein separations: study of size-exclusion media with narrow particle size distributions. J Chromatogr 589:71-85. [184] Bird RB, Stewart WE, Lightfoot EN. 2002. Transport Phenomena. New York: John Wiley and Sons. p. 190. [185] Cao X, Eisenthal R, Hubble J. 2002. Detachment strategies for affinity adsorbed cells. Enzyme Microb Technol 31:153-160. [186] Dismer F, Hubbuch J. 2007. A novel approach to characterize the binding orientation of lysozyme on ion-exchange resins. J Chromatogr A 1149:312-320. [187] Dismer F, Petzold M, Hubbuch J. 2008. Effects of ionic strength and mobile phase pH on the binding orientation of lysozyme on different ion-exchange adsorbents. J Chromatogr A. 1194:11-21 [188] Voitl A, Butte A, Morbidelli M. 2010. Behavior of human serum albumin on strong cation exchange resins: 1. Experimental analysis. J Chromatogr A 1217:5484-5491. [189] Voitl A, Butte A, Morbidelli M. 2010. Behavior of human serum albumin on strong cation exchange resins: II. Model analysis. J Chromatogr A 1217:54925500. [190] Yuan QS, Rosenfield A, Root TW, Klingenberg DJ, Lightfoot EN. 1999. Flow distribution in chromatographic columns. J Chromatogr A 831:149-165. [191] Miyabe K, Guiochon G. 1999. Influence of column radial heterogeneity on peak fronting in linear chromatography. J Chromatogr A 857:69-87. [192] Nakanishi K, Sakiyama T, Imamura K. 2001. On the adsorption of proteins on solid surfaces, a common but very complicated phenomenon. J Biosci Bioeng 91:233-244. [193] Schuck P. 1997. Use of surface plasmon resonance to probe the equilibrium and dynamic aspects of interactions between biological macromolecules. Annu Rev Biophys Biomol Struct 26:541-566. [194] Krägel J, Wüstneck R, Husband F, Wilde PJ, Makievski AV, Grigoriev DO, Li JB. 1999. Properties of mixed protein/surfactant adsorption layers. Colloids Surf, B 12:399-407.  133  [195] Vermeer AWP, Norde W. 2000. CD spectroscopy of proteins adsorbed at flat hydrophilic quartz and hydrophobic Teflon surfaces. J Colloid Interface Sci 225:394-397. [196] Sane SU, Cramer SM, Przybycien TM. 1999. Protein structure perturbations on chromatographic surfaces. J Chromatogr A 849:149-159. [197] Brandes N, Welzl PB, Werner C, Kroh LW. 2006. Adsorption-induced conformational changes of proteins onto ceramic particles: differential scanning calorimetry and FTIR analysis. J Colloid Interface Sci 299:56-69. [198] van der Vegt W, van der Mei HC, Busscher HJ. 1993. Interfacial free energies in protein solution droplets on FEP-teflon by axisymmetric drop shape analysis by profile – IgG versus BSA. J Colloid Interface Sci 156:129-136. [199] Butt H-J. 1996. A sensitive method to measure changes in the surface stress of solids. J Colloid Interface Sci 180:251-260. [200] Maggon K. 2007. Monoclonal antibody “gold rush.” Curr Med Chem 14:19781987. [201]  Shukla AA, Hubbard B, Tressel T, Guhan S, Low D. 2007. Downstream processing of monoclonal antibodies – application of platform approaches. J Chromatogr B 848:28-39.  [202] Lawson K. 2008. Monoclonal therapeutics and companion diagnostic products. BCC Research. Web. 23 June 2011. <http://www.bccresearch.com/report/monoclonal-therapeutics-productsbio016g.html> [203] Farid SS. 2007. Process economics of industrial monoclonal antibody manufacture. J. Chromatogr. B 848:8-18. [204] Boi C. 2007. Membrane adsorbers as purification tools for monoclonal antibody purification. J. Chromatogr. B 848:19-27. [205] Roque ACA, Lowe CR, Taipa MA. 2004. Antibodies and genetically engineered related molecules: production and purification. Biotechnol. Prog. 20:639-654. [206] Roque ACA, Silva CSO, Taipa MA. 2007. Affinity-based methodologies and ligands for antibody purification: advances and perspectives. J. Chromatogr. A 1160:44-55. [207] Zhou JX, Tressel T. 2006. Basic concepts in Q membrane chromatography for large-scale antibody production. Biotechnol Prog 22:341-349.  134  [208] Ochoa-Tapia JA, Whitaker S. 1995. Momentum transfer at the boundary between porous medium and a homogeneous fluid II: comparison with experiment. Int J Heat Mass Transfer 36:2647-2655 [209] Wiedemann K, Roethe A, Radeke K-H, Gelbin D. 1978. The modeling of adsorption-desorption breakthrough cruves using statistical moments. Chem Eng J 16:19-26  135  Appendices Appendix A: Description and program code for the determination of shear profiles from Fluent output There are two components of shear at the membrane surface:     v    z   (A-1)     v r    z   (A-2)   z       zr      The total shear can then be calculated from these components as:        total   z2   zr2  (A-3)  The following MATLAB program was written in order to read velocity vectors calculated by Fluent and determine shear stress profiles at the membrane surface. The program reads the velocities, takes the four data points adjacent to the membrane and fits them to a third order polynomial using the Levenberg-Marquardt method (as seen in the MATLAB function “lev2” shown below) [162].  The derivative of this polynomial is then  determined so as to calculate shear stress.  136  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % This program takes the Fluent calculated velocities from a plane % underneath the rotor, sorts the useful points and fits a 3rd order % polynomial to the four points closest to the membrane using the % Levenberg-Marquardt method (the functions “lev2” and “f12”) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clc clear data = load('alter-P2-60.txt');  % text file must be in proper format % to be read correctly  [d e] = size(data); for i = (1:d) y(i) = data(i,2); z(i) = data(i,3); xvel(i) = data(i,4); yvel(i) = data(i,5); zvel(i) = data(i,6);  % assigns proper names to % velocity data  end j = 1; for i = (2:d) newj = 1; if y(i) == -0.0200 if i < 4 error('Insufficient data points for derivative'); elseif y(i) ~= y(i-1) for k = 1:4 if z(i+1-k) == z(i-k) r(k,j) = z(i+1-k)*1000; % collates only the tanvel(k,j) = xvel(i+1-k)*1000; % 4 data points radvel(k,j) = zvel(i+1-k)*1000; % above the axvel(k,j) = yvel(i+1-k)*1000; % membrane for ht(k,j) = y(i+1-k)*1000+20; % derivative else % calculation newj = 2; % also must end % translate to end % appropriate form if newj == 1 % i.e 0 = membrane j = j+1; % surface end end end end [d2 e2] = size(ht); x0 = [10 1000 10000 0];  % initial guesses for polynomial fitting  137  for i = 1:e2 dtandy = lev2(ht(:,i),tanvel(:,i),x0); % lev2 fits 4 data pts to a draddy = lev2(ht(:,i),radvel(:,i),x0); % 3rd order polynomial daxdy = lev2(ht(:,i),axvel(:,i),x0); shears(i) = 0.001*(dtandy(2)^2 + draddy(2)^2)^0.5; end radpos = r(1,:)/1000; outaxvel = axvel(1,:)/1000;  % also collates fluid velocity through % membrane to ensure it is constant  save shear.out radpos shears -ASCII save axvel.out radpos outaxvel -ASCII % output is in m for radpos Pa for shears and m/s for axvel %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “lev2” fits data to a 3rd order polynomial using the % Levenberg-Marquardt method. This function calls another function, % “f12”. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function A = lev2(x,y,a) lambda = 1000; dela = a; [m n] = size(a); xi = f12(x,y,a); while (abs(max(dela))>1e-8 || abs(min(dela))>1e-8) for k = 1:n apos(k) = a(k) + 1e-3; aneg(k) = a(k) - 1e-3; for kk = 1:n if kk == k dumpos(kk) dumneg(kk) else dumpos(kk) dumneg(kk) end end  = apos(kk); = aneg(kk); = a(kk); = a(kk);  beta(k) = -1/2*(f12(x,y,dumpos) - f12(x,y,dumneg))/2e-3; end for j = 1:n for p = 1:n  138  if p == j for q = 1:n if q ~= p dumpos2(q) dumneg2(q) else dumpos2(q) dumneg2(q) end end  = a(q); = a(q); = apos(p); = aneg(p);  alpha(j,p) = 1/2*(f12(x,y,dumpos2) + f12(x,y,dumneg2) -… 2*f12(x,y,a))/1e-6; else for r = 1:n if r == p dumpos3(r) = apos(r); dumneg3(r) = aneg(r); dumposneg(r) = aneg(r); dumnegpos(r) = apos(r); elseif r == j dumpos3(r) = apos(r); dumneg3(r) = aneg(r); dumposneg(r) = apos(r); dumnegpos(r) = aneg(r); else dumpos3(r) = a(r); dumneg3(r) = a(r); dumposneg(r) = a(r); dumnegpos(r) = a(r); end end alpha(j,p) = 1/2*(1/(4e-6)*(f12(x,y,dumpos3)-… f12(x,y,dumnegpos)- f12(x,y,dumposneg) +… f12(x,y,dumneg3))); end end end aprime = alpha + lambda*eye(n); dela = aprime\beta'; a2 = a + dela'; xinew = f12(x,y,a2); if xinew < xi a = a2; xi = xinew; lambda = lambda/5; else lambda = lambda*5; end end A = a;  139  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “f12” defines the 3rd order polynomial to which the data % is fit and calculates the chi squared values for the % Levenberg-Marquardt method. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function sum = f12(u,v,a) sum = 0; [d e] = size(u); for i = 1:d f12 = (v(i) - (a(1) + a(2)*u(i) + a(3)*u(i)^2 + a(4)*u(i)^3))^2; newsum = sum + f12; sum = newsum; end  140  Appendix B: Description of the similarity solution for a two-domain model in the absence of suction through the bottom wall The physical configuration of the simplified CSAF system for which a two-domain model could be derived and an analytical solution obtained for tangential fluid velocities within the rotor chamber (fluid domain D1) and membrane stack (porous domain D2) is shown in Figure B.1. The equations of motion are defined in Section 1.4 for both the fluid domain (Equations 1-1 and 1-2) and the porous domain (Equations 1-11 and 1-12). The momentum equations can be expanded in cylindrical coordinates to yield: r-momentum equations:  D1:  m   u r u v r u 2 u  ur    uz r    r r  r z   (B-1) 2 2   1  u  u  u    p ru r   12 2r  22   2r   m u r  m    r K r  z   r   r  r r  z    D2 D1  h(r)   Ho  r  Um Figure B.1  Physical configuration of the 2-domain similarity solution for CSAF hydrodyamics. A conical rotor with an angle of incidence  and a localized height, h(r), above the porous medium of thickness Ho, is rotated at a frequency . The axial velocity at z = 0 is specified as Um.  141      v r D2:    2 v r v v r v v     v z r   r r  r z  2 2   1   p rvr   12  v2r  22 v   v2r      r r  z   r   r  r r  (B-2)  momentum equations:  D1:  u    u u u u r u    uz     ur  m  r r  r z    1   2u  2u ru   12 2  22 u r  2   r  z  r   r  r r      vr   D2:   m u  K   v v v v r v v     vz    r r  r z  2   1  2 v r  2 v   1  v rv   2 2  2     r  z 2   r   r  r r  (B-3)  (B-4)  z-momentum equations:  D1:  u    u z u u z   uz z    ur  m  r r  z   1   u z  1  2 u z  2 u z   m  p  m    uz r  2  2 z K z 2   r r  r  r       vr D2:    v z v v z v    vz z   r r  z   1   v z  1  2 v z  2 v z  p     2  r  2 2 z z   r r  r  r   (B-5)  (B-6)  The boundary conditions are: no slip at the rotating wall and lower boundary of the porous medium; equivalent fluid velocities at the interface of the two domains and a discontinuity in the shear stress at the interface which, as proposed by Ochoa-Tapia and Whitaker [206], is inversely proportional to the permeability of the porous medium: at the rotating wall  vr = vz = 0 v= r  (B-7) 142  at the membrane  at the lower boundary  u(r,,z) = v(r,,z)  (B-8)  P+ = P-  (B-9)  1 u r v r Bu r    m z z K  (B-10)  Bu 1 u v     m z z K  (B-11)  ur = u= 0uz = Um  where B is an empirical constant determined, in this case, to be 0.15 through comparison to CFD results. Using a similarity solution the equations of motion can be reduced to a more solvable form. Following the procedure outlined by von Karman [130], streamlines in domains D1 and D2 are postulated to take the form:    z   H ( r )    (B-12)   z    H r    (B-13)   1 r , z    U m r 2 g1  1 2   2 r , z    U m r 2 g 2  1 2  The fluid velocities in both domains can then be written as functions of these streamline functions as well as dimensionless tangential velocity functions 1 and 2:  ur   1  1 r , z  1 U m r    z    g1  r z 2 H r  z   H r    z   u θ  r1   H r    (B-14)  (B-15)  143   d    z     r 2 H r   g1    z  1 dr r   H r    1  1 r , z    uz    U m  g 1   2 r r H r    H r   2       vr   1  2 r , z  1    z     U m r  g 2  r z 2 z   H r    (B-17)   z   vθ  r2   H r    (B-18)   d    z     r 2 H r   g 2    z  1 dr r   H r    1  2 r , z    vz    U m  g 2   2   r r H r H r    2       where H(r) = Ho + h(r).  (B-16)  (B-19)  These assumed functions satisfy both the integral and  differential continuity equations. After substitution of the above velocity expressions into Equations B-1 through B-6, elimination of the pressure terms through cross differentiation of the r- and z-momentum equations in each domain, non-dimensionalization and omission of terms of order H r  the following set of four ordinary differential equations is obtained: 2 ε H ξ  d  d  d 2 αH ξ    λ1 η g 1 η  g 1 η λ1 η   2 λ1 η  m λ1 η  0  ε m Rrot  K  dη   dη   dη  (B-20)  d  d  d 2 αH ξ    λ2 η g 2 η  g 2 η λ2 η   2 λ2 η  0  Rrot   dη   dη   dη  (B-21)   d4  4β 2 H ξ 3  d3  ε H ξ  d  αH ξ   4 g 1 η   3 g 1 η  m         λ η λ η  g η 1 1 1   3 K  dη  ε m Rrot  dη  αε m Rrot  dη   2   d2   2 g 1 η  0  dη   (B-22)  144   d4  4β 2 H ξ 3  d3  d  αH ξ   4 g 2 η   3 g 2 η  0         λ η λ η  g η 2 2 2 3   αR rot Rrot  dη   dη   dη   (B-23)  The non-dimensional coordinates  and  are given by: ξ  r Rrot  η  z  (B-24)  (B-25)  H r   and  and  are the Reynolds numbers for flow through the membrane and within the rotor chamber, respectively, defined by Equations 2-9 and 2-10. Domain 1 is defined as: (,)=[0,1],[0,] and Domain 2 as: (,)=[0,1],[1] where is the ratio Ho/H and represents the interface between the fluid and porous domains. This formulation is valid in the case of nearly constant rotor height over small radial distances. In the limit where there is no suction through the bottom wall (0)Equations B-20 and B-21 reduce to:  ε m H   d2   λ η  λ1 η  0 1 K dη 2  (B-26)  d2 λ2 η   0 dη 2  (B-27)  2  The boundary conditions are then given by:  λ1 0   0 λ2 1  1  λ1 ζ   λ2 ζ   K   λ1 ζ   ε m λ2 ζ   ε m B 2    H     (B-28) 1/2  λ1 ζ   145  Equations B-26 and B-27, coupled with the boundary conditions, can be solved analytically at a specified value of  to yield the following expressions for the tangential velocities in the porous and fluid domains:  1    A1e   m     Da      A2 e      m   Da    (B-29)  2    A3  A4  (B-30)  where,  A1   A2   A3     m Dae  2      m B1     Da 1  e        m Da     m   Da         2             1   1  e m           m Dae     m   Da   m Da              2 m   2 m             Da Da       m B1     Da 1  e    m 1   1  e             m     m e   2 m     Da     B m   m   2 m   2 m             Da Da       m B1     Da 1  e    m 1   1  e              2 m   2 m            Da   Da      m  1  e   B   Da 1  e  m           A4   2 m   2 m            Da   Da      m B1     Da 1  e    m 1   1  e               (B-32)    B   (B-31)     (B-33)    (B-34)    146  Appendix C: Description and program code for the one-domain similarity solution The one-domain similarity solution strategy is described in Section 2-3 and is derived in a similar fashion as the two-domain solution (see Appendix B). Equations 2-4 and 2-5 have the potential to become stiff and so for ease of solution they can be further scaled to remove the Reynolds number terms and  can be stretched:    g    b 2 ˆ  (C-1)   bgˆ  (C-2)      (C-3)  b  Substituting these expressions into Equations 2-4 and 2-5 yields:   d ˆ  d2 ˆ h   ˆ  d     ˆ ˆ           g   g       0 2  Rrot  d  d   d   (C-4)   d4  h   d3  4h 3 ˆ  d ˆ   4 gˆ    gˆ   3 gˆ           0 3 Rrot  d   d  Rrot  d   (C-5)  with the boundary conditions: gˆ b   gˆ b   0 gˆ 0  ˆ 0  0  ˆb     b2  gˆ 0      (C-6)  b  The set of equations can be written in the following form for easy translation to a RungeKutta numerical solution [160]:  147  ˆ    ˆ1  ˆ  g   gˆ 1     gˆ 2   gˆ 3         ˆ1     h  ˆgˆ  gˆˆ  1 1   Rrot   gˆ 1     gˆ 2   gˆ 3     3  4h  ˆˆ  h  gˆgˆ  1 3  R 3  Rrot      (C-7)  where,  ˆ   ˆ1   gˆ   gˆ 1    2 gˆ   gˆ 2   2  3 gˆ   gˆ 3   3  (C-8)  Furthermore, the two components of shear at the membrane surface can be written in terms of the functions    and g   as:    z     r    h r        hr       (C-9)    z    U m rg     h r 1     zr      2  hr       (C-10)     z  After all variable transformations, non-dimensionalization and scaling these shear equations take the form:  148     rˆ1  b 3    hr     z      (C-11)   1 U m rgˆ 3  b3   zr     2  2 hr    (C-12)    The MATLAB code for the one-domain similarity solution is shown below.  149  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % This program solves the one-domain numerical problem using the built % in ode solver ode15s based on scaled and stretched variables. % The program calls the function “a1eqn2” which contains the equations % to be solved. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clc clear global R h delv2 delv5 delv6 eps =  = 1e-4; = 1e-4; = 1e-4; 1e-3;  % tolerance of 0.1% between upper boundary % conditions and calculated values  nmax = 10000;  % maximum number of iterations allowed  R = 35e-3; dens = 998.2; mu = 1.00e-3;  % Rotor radius % density of water % viscosity of water  omega = 1.047;  % rotational velocity of rotor in rad/sec  Vw = 3.65e-5; r = 5e-3;  % fluid flow through membrane % radial position  xi = r/R;  % non-dimensional radial position  h = r*0.06993 + 0.0002; % height at radial position, in this case the % rotor is a cone of angle 4 degrees with a gap % of 0.0002 m alpha = dens*Vw*R/mu; beta = dens*omega*R^2/mu;  % Reynolds number alpha % Reynolds number beta  b = 5;  % stretching factor  vbot = [0;78;-alpha/b;0;10;-9]; % initial guesses for velocity vectors vtop = [beta/b^2;1;0;0;-1.9;14]; % at membrane and rotor surfaces F = [2;2;2];  % difference between upper BCs and solution  count = 0; while max(abs(F)) > (eps)  % loop will run until solution % matches upper BCs [eta,v1] = ode15s(@a1eqn2,[0 c],vbot); % solves the system of [s,q] = size(v1); % odes using built in % solver  150  F(1) = vtop(1) - v1(s,1); F(2) = vtop(3) - v1(s,3); F(3) = vtop(4) - v1(s,4);  % determines discrepancy between % upper BC and solution  dumvbot2 = vbot + [0;delv2*vbot(2);0;0;0;0]; dumvbot5 = vbot + [0;0;0;0;delv5*vbot(5);0]; dumvbot6 = vbot + [0;0;0;0;0;delv6*vbot(6)];  % % % % %  writes new bottom vectors with slightly different values  % solves with new bottom vectors [eta2,dumv2] = ode15s(@a1eqn2,[0 b],dumvbot2); [s2,q2] = size(dumv2); [eta5,dumv5] = ode15s(@a1eqn2,[0 b],dumvbot5); [s5,q5] = size(dumv5); [eta6,dumv6] = ode15s(@a1eqn2,[0 b],dumvbot6); [s6,q6] = size(dumv6); % use the solutions found above to determine the differential % of the F vector with respect to the bottom velocities dFdv(1,1) = 1/(delv2*vbot(2))*(vtop(1) - dumv2(s2,1) - F(1)); dFdv(1,2) = 1/(delv5*vbot(5))*(vtop(1) - dumv5(s5,1) - F(1)); dFdv(1,3) = 1/(delv6*vbot(6))*(vtop(1) - dumv6(s6,1) - F(1)); dFdv(2,1) = 1/(delv2*vbot(2))*(vtop(3) - dumv2(s2,3) - F(2)); dFdv(2,2) = 1/(delv5*vbot(5))*(vtop(3) - dumv5(s5,3) - F(2)); dFdv(2,3) = 1/(delv6*vbot(6))*(vtop(3) - dumv6(s6,3) - F(2)); dFdv(3,1) = 1/(delv2*vbot(2))*(vtop(4) - dumv2(s2,4) - F(3)); dFdv(3,2) = 1/(delv5*vbot(5))*(vtop(4) - dumv5(s5,4) - F(3)); dFdv(3,3) = 1/(delv6*vbot(6))*(vtop(4) - dumv6(s6,4) - F(3)); del=dFdv\-F; vbot(2) = vbot(2) + del(1); vbot(5) = vbot(5) + del(2); vbot(6) = vbot(6) + del(3);  % write new velocity vector at % membrane surface  count = count + 1; if count > nmax error('solution not converging'); end  % ensures that an % infinite loop has not % been created  end % with the determined velocity vectors return non-stretched and non% scaled values for lambda, g and g prime lambda = b^2*v1(:,1)/beta; g = b*v1(:,3)/alpha; gp = b^2*v1(:,4)/alpha; tzt = -mu*r*omega*v1(1,2)*b^3/(beta*h); trz = mu*r*Vw/2*v1(1,5)/alpha*b^3/(h)^2;  % calculates shear components  151  shear = (tzt^2 + trz^2)^0.5;  % calculates shear  eta = eta/c; eta = eta'; lambda = lambda'; g = g'; gp = gp';  % ensures the vectors are in the correct % orientation  save veldata.out eta lambda g gp -ASCII  % saves output  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “a1eqn2” is used in the above program; it defines the % set of equations to be solved % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dv = a1eqn2(eta,y) global R h dv = zeros(6,1); dv(1) = y(2); dv(2) = -h/R*(y(1)*y(4) - y(3)*y(2)); dv(3) = y(4); dv(4) = y(5); dv(5) = y(6); dv(6) = 4*h^3/R^3*y(1)*y(2) + h/R*y(3)*y(6);  152  Appendix D: Description of modeling equations for a CSTR with multiple inflows The sum of volumetric flow rates Qi through the individual zones of the ZRM equals the overall volumetric flow rate Q through the membrane stack. Flow fractions i = Qi/Q can therefore be defined which are independent of the overall volumetric flow rate Q and which contain one degree of freedom less than the corresponding volumetric flow rates Qi. Together, the flow fractions and the interstitial velocities are required to determine the size of each membrane zone since a small zone with fast interstitial flow can have the same volumetric flow rate as a large zone with slow interstitial flow. In a network of CSTRs the number of outflows from a given tank is irrelevant to the mathematical description.  If one tank feeds another in the physical model, the  corresponding equations are coupled by taking the concentration within the feeding tank as the inlet concentration of the fed tank. The same rule holds when a tank feeds a membrane zone. Downstream of the membrane stack however, the model includes tanks with multiple inflows (see Figure 4.1). A tank with j inflows is mathematically described by an instance of Equation D-1: m c m c CSTR 1 in , j   c CSTR  t j 1  j j 1  j CSTR  (D-1)  Here cCSTR denotes the concentration in the respective tank, c CSTR the corresponding inlet in concentration and  the residence time. A set of flow fractions must be defined for each tank with more than one inlet. Figure D.1 illustrates this situation in more detail and represents the collector tanks from Figure 4.1a. The set 1 to 3 with 1 + 2 + 3 = 1 refers to the flow fraction through each zone, while the sets 1b,1 and 1b,2 with 1b,1 + 1b,2 = 1 as well as 2b,1 and 2b,2 with  2b,1 + 2b,2 = 1 refer to the inlets of tanks 1b and 2b, respectively. An instance of Equation D-1 is set up for each of these tanks with residence times 1/1b,1 + 1/1b,2 = 1/1b  153  3b  2b,1 2b 1b  2b,2  1b,2  1b,1  Figure D.1  Junction of flows through three zones in a cascade of CSTRs.  and 1/2b,1 + 1/2b,2 = 1/2b, and with flow fractions 1b,i = 1b/1b,i and 2b,i = 2b/2b, for 1  i  2 . The flow fractions through the membrane zones are then uniquely determined  by 1 = 1b,1, 2 = 1b,2 · 2b,1 and 3 = 1b,2 · 2b,2. Analogous relations can be derived for complex networks with more membrane zones or with exchange flows between the tanks.  154  Appendix E: Program code for the zonal rate model The following MATLAB code was written by Dr. Eric von Lieres and represents the programming framework for the ZRM. The solution of the set of differential-algebraic equations that constitute the model relies on the built-in ode15s solver. The efficiency of this solver is greatly increased by using custom functions to re-shape the system Jacobian. These custom functions are not included below for proprietary reasons and so while the following code can be used to optimize parameter values for the ZRM it does so orders of magnitude slower than the times listed in Section 4.2.2.  155  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “parameter” defines the required model parameters. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function parameter global te ts ze zs rt at u Ez ec ep sv Kf c0 t0 dbl dbr wo mo fd n it qs zm zr col_area zm = 21; switch case case case case case case case case case case case case case case case case case end  zm 10; 11; 13; 14; 15; 16; 21; 22; 23; 24; 25; 26; 31; 33; 34; 36; 44;  % defines what model configuration to implement zr zr zr zr zr zr zr zr zr zr zr zr zr zr zr zr zr  = = = = = = = = = = = = = = = = =  1; 100; [175 25]; [1 2 3]; [1 2 3 4]; [1 2 3 4 5 6]; [139.69 79.37 278.01]; [1 2 3 4]; [1 2 3 4 5]; [1 2 3 4 5]; [1 2 3 4 5 6 7]; [1 2 3 4 5]; [1 2 3 4 5]; [1 2 3 4 5 6 7 8]; [1 2 3 4 5 6 7]; [1 2 3 4 5 6 7 8]; [1 2 3 4 5 6 7 8 9];  %tank network residence %times  % domain ts = 200; zs = 200;  % time steps (for output) % spatial knots  % Danckwerts boundary condition (0 no/ 1 yes) - integer values >1 % specify order of corresponding (non-central) differential quotient dbl = 1; dbr = 1; % solver rt = 1e-4; at = 1e-9;  % relative tolerance % absolute tolerance  n = 1; fd = 0; it = 'scl';  % Number of components % Film-Diffusion (0 no/ 1 yes) % Isotherm ('sma' Steric Mass Action/ 'scl'  156  % Single Component Langmuir/ 'mcl' Multi % Component Langmuir/'bil' Bi-Langmuir % /'spr' Spreading/'sprS' Simplified Spreading % times te = 1200; t0 = te;  % simulation time  % geometry ze = 2.2e-3; ec = 0.70; col_area = 2.206e-3;  % column length (m) % column porosity % frontal column area (m^2)  % transport % u is interstitial (m/s)  first number is flowrate in mL/min  u = 1.5 * (1/60) * (1/100^3) * (1/col_area) * (1/ec); Ez = 7e-11; %axial dispersion coeff. m^2/s % % % % % % % %  Film mass transfer resistance epapp is the chosen percentage of void volume that is stagnant liquid film ecapp is the given column voidage ep translates apparent voidages to what the model must see dp is the pore radius sv is the surface to volume ratio of the pore Kf is the product of the film diffusion coeff and sv  if fd > 0 epapp = 0.01; ecapp = 0.7; ep = ecapp*epapp/(1 – ecapp + epapp*ecapp); dp = 0.8e-6; sv = 4/dp; Kf = 1 * sv; u = 1.5 * 7.55515e-6/ec; end % Imposed concentrations c0 = 1; % Isotherm parameters global ni qstar qm qm2 ka kd beta ka2 kd2 k12 k21 s v d ka1 kd1 ni = 1; if strcmp(it,'scl') qm = 0; ka = 0; kd = 0; elseif strcmp(it,'mcl') n = 2; ni = 2; c0 = [1 1];  157  qm = [0 0]; ka = [0 0]; kd = [0 0]; Ez = [Ez 0]; elseif strcmp(it,'bil') qm = 150; qm2 = 150; ka = 9.8e-3; kd = 2.6e-5; ka2 = 1e-3; kd2 = 5e-4; elseif strcmp(it,'sma') n = 2; ni = 2; c0 = [0.875 1]; ka1 = [0 0]; kd1 = [0 0]; s = [0 0]; v = [0 0]; d = 0; Ez = [Ez Ez]; elseif strcmp(it,'spr') qm = 0; ka = 0; kd = 0; ka2 = 0; kd2 = 0; k12 = 0; k21 = 0; beta = 0; elseif strcmp(it,'sprS') qstar = 400; ka = 1.96e-2; kd = 3.5e-4; beta = 2.43; k12 = 0; k21 = 0; end  158  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “optimize_zonal” is used to optimize parameters. It % calls the function “residual_zonal” which simulates a breakthrough % curve from given parameter values then calculates the residual % between this result and given BTC data. “optimize_zonal” minimizes % this residual using the built in function lsqonlin and then returns % the optimized parameter values. % In this version of the function the parameters that are being % optimized are the isotherm parameters ‘ka’ ‘kd’ and ‘qm’ for four % flowrates: 1.5, 5, 10 and 20 mL/min. ‘ka’ and ‘kd’ are flowrate % independent whereas ‘qm’ varies with flow. Only experimental data % under binding conditions is being used for parameter optimization. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function optimize_zonal clear global; global data c0 zm nc sv ec global nd  % necessary % defined  % prerequisites parameter;  % sets the parameter values from “parameter” function  discrete;  % calls the “discrete” function seen below  % load data data(1).file data(1).comp  = 'avg-binding-1.5mLmin-1mgmL.mat'; = 1; % 0 = uses computation from previous dataset, % 1 = performs computation  data(2).file data(2).comp  = 'avg-binding-5mLmin-1mgmL.mat'; = 1;  data(3).file data(3).comp  = 'avg-binding-10mLmin-1mgmL.mat'; = 1;  data(4).file data(4).comp  = 'avg-binding-20mLmin-1mgmL.mat'; = 1;  nd = length(data); for i = 1:nd load(data(i).file); data(i).tm = tm; data(i).Cm = Cm; end  159  % initial guesses for the parameters to be optimized param = [0.005 0.005 300 300 300 300]; param_log = log10(param);  % translates param into logarithmic scale  sum(residual_zonal(param_log).^2)  % determines the initial residual  % optimized values for the parameters are determined by using the built % in function lsqnonlin which minimizes the residual calculated by the % function “residual_zonal”. param_opt_log = lsqnonlin(@residual_zonal,param_log,0*param_log-… 16,0*param_log+16,optimset('display','iter','DiffMinChange',1e-… 4,'MaxFunEvals',Inf,'TolX',1e-4),out); dev_opt res_lsq res_abs res_max  = = = =  residual_zonal(param_opt_log); sum(dev_opt.^2); sum(abs(dev_opt)); 100*max(abs(dev_opt));  % sum of square of residuals % sum of absolute residuals % maximum deviation from data  param_opt = 10.^param_opt_log;  % the optimized values % of the parameters  param_opt = 10.^param_opt_log;  % translates out of log  param_opt_log param_opt res_lsq  % outputs optimized parameter % values and residual  160  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “discrete” is called by “optimize_zonal” and % “residual_zonal” in order to define the matrices required for % solution. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function discrete global te ze ra ts zs rs nc zm n zr Ez fd Kf wc ep dbl it % necessary global t z r dt dz dr zi nt % defined %nc is the number of zones %nt is the number of tanks %zi is the identity of the tanks that inlet to the membrane stack switch case case case case case case case case case case case case case case case case case end  zm 10; 11; 13; 14; 15; 16; 21; 22; 23; 24; 25; 26; 31; 33; 34; 36; 44;  nc nc nc nc nc nc nc nc nc nc nc nc nc nc nc nc nc  = = = = = = = = = = = = = = = = =  1; 1; 1; 1; 1; 1; 2; 2; 2; 2; 2; 2; 3; 3; 3; 3; 4;  nt nt nt nt nt nt nt nt nt nt nt nt nt nt nt nt nt  = = = = = = = = = = = = = = = = =  1; 2; 2; 4; 4; 6; 4; 4; 4; 6; 6; 4; 6; 6; 8; 6; 10;  zi zi zi zi zi zi zi zi zi zi zi zi zi zi zi zi zi  = = = = = = = = = = = = = = = = =  1; 1; 1; 2; 2; 3; [1 [1 [1 [2 [2 [1 [1 [1 [2 [1 [2  2]; 2]; 2]; 3]; 3]; 2]; 2 3]; 2 3]; 3 4]; 2 3]; 3 4 5];  % Time discretization dt = te/(ts - 1); t = 0:dt:te; % Column discretization dz = ze/(zs - 1); z = 0:dz:ze; global global global global global global  dz dbl dbr; M2; d0 d1 d2; A0 A1 A2; B0 B1 B2; C0 C1 C2;  % necessary % defined  161  % Derivatives for Danckwerts boundary condition A0 = spalloc(zs,zs,zs); switch case case case end  dbl 1; A0(1,1:2) = [-1 1]/dz; 2; A0(1,1:3) = [-3 4 -1]/(2*dz); 3; A0(1,1:4) = [-11 18 -9 2]/(6*dz);  switch case case case end  dbr 1; A0(zs,zs-1:zs) = [-1 1]/dz; 2; A0(zs,zs-2:zs) = [1 -4 3]/(2*dz); 3; A0(zs,zs-3:zs) = [-2 9 -18 11]/(6*dz);  % First derivative A0(2:zs-1,1:zs-2) = A0(2:zs-1,1:zs-2) - 1/dz*speye(zs-2); A0(2:zs-1,2:zs-1) = A0(2:zs-1,2:zs-1) + 1/dz*speye(zs-2); if not(dbr); A0(zs,zs-1:zs) = [-1 1]/dz; end % Second derivative M2 = spalloc(zs,zs,(2*mo+1)*zs); M2(1,1:3) = [1 -2 1]/dz^2; for i = 2:zs-1 M2(i,i-1:i+1) = [1 -2 1]/dz^2; end M2(zs,zs-2:zs) = [1 -2 1]/dz^2; global n wc Ez fd Kf zm zr nc nt;  % necessary  EZ = repmat(Ez,[zs,1,nc]); if fd > 0 KF = repmat(Kf,[zs,1,nc]); end % Jacobian for transport equations global JFC JFF % defined JFC = spalloc(0,0,0); JFF = spalloc(0,0,0); if fd > for i JFC JFF end end  0 = 1:n = blockdiag( JFC , 1/ep*Kf(i)*eye(zs) ); = blockdiag( JFF , -1/ep*Kf(i)*eye(zs) );  162  % Boundaries if not(dbl) for i = 1:n j = (i - 1)*zs + 1; JFC(j,:) = 0; end end % Zones global JTT JTC; % defined JTT = spalloc(nt,nt,0); JTC = spalloc(n*nt,nc*n*zs,0); switch zm case 10; JTT(1,1) = case 11; JTT(1,1) = JTT(2,2) = case 13; JTT(1,1) = JTT(2,2) = case 14; JTT(1,1) = JTT(2,1) = JTT(2,2) = JTT(3,3) = JTT(4,3) = JTT(4,4) = case 15; JTT(1,1) = JTT(2,1) = JTT(2,2) = JTT(3,3) = JTT(4,3) = JTT(4,4) = case 16; JTT(1,1) = JTT(2,1) = JTT(2,2) = JTT(3,2) = JTT(3,3) = JTT(4,4) = JTT(5,4) = JTT(5,5) = JTT(6,5) = JTT(6,6) = case 21; JTT(1,1) = JTT(2,1) = JTT(2,2) = JTT(3,3) = JTT(4,3) = JTT(4,4) = case 22; JTT(1,1) = JTT(2,1) = JTT(2,2) = JTT(3,3) = JTT(4,3) = JTT(4,4) = case 23; JTT(1,1) =  -1/zr(1); -1/zr(1); -1/zr(1); -1/zr(1); -1/zr(2); -1/zr(1); 1/zr(2); -1/zr(2); -1/zr(2); 1/zr(3); -1/zr(3); -1/zr(1); 1/zr(2); -1/zr(2); -1/zr(3); 1/zr(4); -1/zr(4); -1/zr(1); 1/zr(2); -1/zr(2); 1/zr(3); -1/zr(3); -1/zr(4); 1/zr(5); -1/zr(5); 1/zr(6); -1/zr(6); -1/zr(2)-1/zr(3); 1/zr(1); -1/zr(1); -1/zr(1); 1/zr(3); -1/zr(2)-1/zr(3); -1/zr(1); 1/zr(2); -1/zr(2); -1/zr(2); 1/zr(4); -1/zr(3)-1/zr(4); -1/zr(1);  163  case 24;  case 25;  case 26;  case 31;  case 33;  case 34;  JTT(2,1) JTT(2,2) JTT(3,3) JTT(4,3) JTT(4,4) JTT(1,1) JTT(2,1) JTT(2,2) JTT(3,1) JTT(3,3) JTT(4,4) JTT(5,5) JTT(6,4) JTT(6,5) JTT(6,6) JTT(1,1) JTT(2,1) JTT(2,2) JTT(3,1) JTT(3,3) JTT(4,4) JTT(5,5) JTT(6,4) JTT(6,5) JTT(6,6) JTT(1,1) JTT(2,1) JTT(2,2) JTT(3,3) JTT(4,3) JTT(4,4) JTT(1,1) JTT(2,1) JTT(2,2) JTT(3,2) JTT(3,3) JTT(4,4) JTT(5,4) JTT(5,5) JTT(6,5) JTT(6,6) JTT(1,1) JTT(2,1) JTT(2,2) JTT(3,2) JTT(3,3) JTT(4,4) JTT(5,4) JTT(5,5) JTT(6,5) JTT(6,6) JTT(1,1) JTT(2,1) JTT(2,2) JTT(3,1) JTT(3,3) JTT(4,1)  = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =  1/zr(2); -1/zr(2); -1/zr(3); 1/zr(5); -1/zr(4)-1/zr(5); -1/zr(1); 1/zr(2); -1/zr(2); 1/zr(3); -1/zr(3); -1/zr(2); -1/zr(3); 1/zr(4); 1/zr(5); -1/zr(4)-1/zr(5); -1/zr(1); 1/zr(2); -1/zr(2); 1/zr(3); -1/zr(3); -1/zr(4); -1/zr(5); 1/zr(6); 1/zr(7); -1/zr(6)-1/zr(7); -1/zr(1); 1/zr(2); -1/zr(2); -1/zr(3); 1/zr(5); -1/zr(4)-1/zr(5); -1/zr(4)-1/zr(5); 1/zr(2)+1/zr(3); -1/zr(2)-1/zr(3); 1/zr(1); -1/zr(1); -1/zr(1); 1/zr(3); -1/zr(2)-1/zr(3); 1/zr(5); -1/zr(4)-1/zr(5); -1/zr(1); 1/zr(2); -1/zr(2); 1/zr(3); -1/zr(3); -1/zr(4); 1/zr(6); -1/zr(5)-1/zr(6); 1/zr(8); -1/zr(7)-1/zr(8); -1/zr(1); 1/zr(2); -1/zr(2); 1/zr(3); -1/zr(3); 1/zr(4);  164  JTT(4,4) = -1/zr(4); JTT(5,5) = -1/zr(2); JTT(6,6) = -1/zr(3); JTT(7,7) = -1/zr(4); JTT(8,5) = 1/zr(5); JTT(8,6) = 1/zr(6); JTT(8,7) = 1/zr(7); JTT(8,8) = -1/zr(5)-1/zr(6)-1/zr(7); case 36; JTT(1,1) = -1/zr(1); JTT(2,1) = 1/zr(2); JTT(2,2) = -1/zr(2); JTT(3,2) = 1/zr(3); JTT(3,3) = -1/zr(3); JTT(4,4) = -1/zr(4); JTT(5,4) = 1/zr(6); JTT(5,5) = -1/zr(5)-1/zr(6); JTT(6,5) = 1/zr(8); JTT(6,6) = -1/zr(7)-1/zr(8); case 44; JTT(1,1) = -1/zr(1); JTT(2,1) = 1/zr(2); JTT(2,2) = -1/zr(2); JTT(3,1) = 1/zr(3); JTT(3,3) = -1/zr(3); JTT(4,1) = 1/zr(4); JTT(4,4) = -1/zr(4); JTT(5,1) = 1/zr(5); JTT(5,5) = -1/zr(5); JTT(6,6) = -1/zr(2); JTT(7,7) = -1/zr(3); JTT(8,8) = -1/zr(4); JTT(9,9) = -1/zr(5); JTT(10,6) = 1/zr(6); JTT(10,7) = 1/zr(7); JTT(10,8) = 1/zr(8); JTT(10,9) = 1/zr(9); JTT(10,10) = -1/zr(6)-1/zr(7)-1/zr(8)-1/zr(9); end for i = 1:n i1 = (i-1)*nt; i2 = (i-1)*zs; switch case case case case case case  zm 11; 13; 14; 15; 16; 21;  case 22; case 23; case 24; case 25;  JTC(i1+2,i2+ JTC(i1+2,i2+ JTC(i1+3,i2+ JTC(i1+3,i2+ JTC(i1+4,i2+ JTC(i1+3,i2+ JTC(i1+4,i2+ JTC(i1+3,i2+ JTC(i1+4,i2+ JTC(i1+3,i2+ JTC(i1+4,i2+ JTC(i1+4,i2+ JTC(i1+5,i2+ JTC(i1+4,i2+  zs) zs) zs) zs) zs) n*zs+zs) zs) n*zs+zs) zs) n*zs+zs) zs) zs) n*zs+zs) zs)  = = = = = = = = = = = = = =  1/zr(1); 1/zr(2); 1/zr(2); 1/zr(3); 1/zr(4); 1/zr(1); 1/zr(2); 1/zr(2); 1/zr(3); 1/zr(3); 1/zr(4); 1/zr(2); 1/zr(3); 1/zr(4);  165  JTC(i1+5,i2+ n*zs+zs) case 26; JTC(i1+3,i2+ zs) JTC(i1+4,i2+ n*zs+zs) case 31; JTC(i1+4,i2+2*n*zs+zs) JTC(i1+5,i2+ n*zs+zs) JTC(i1+6,i2+ zs) case 33; JTC(i1+4,i2+2*n*zs+zs) JTC(i1+5,i2+ n*zs+zs) JTC(i1+6,i2+ zs) case 34; JTC(i1+5,i2+ zs) JTC(i1+6,i2+ n*zs+zs) JTC(i1+7,i2+2*n*zs+zs) case 36; JTC(i1+4,i2+ zs) JTC(i1+5,i2+ n*zs+zs) JTC(i1+6,i2+2*n*zs+zs) case 44; JTC(i1+6,i2+ zs) JTC(i1+7,i2+ n*zs+zs) JTC(i1+8,i2+2*n*zs+zs) JTC(i1+9,i2+3*n*zs+zs) end end  = = = = = = = = = = = = = = = = = = =  1/zr(5); 1/zr(3); 1/zr(4); 1/zr(1); 1/zr(2); 1/zr(4); 1/zr(4); 1/zr(5); 1/zr(7); 1/zr(2); 1/zr(3); 1/zr(4); 1/zr(4); 1/zr(5); 1/zr(7); 1/zr(2); 1/zr(3); 1/zr(4); 1/zr(5);  % Repetitions H1 = JFC; JFC = []; H2 = JFF; JFF = []; for i = 1:nc JFC = blockdiag(JFC,H1); JFF = blockdiag(JFF,H2); end H1 = JTT; JTT = []; for i = 1:n JTT = blockdiag(JTT,H1); end % Define matrices of the isotherm parameters global ni qm ka1 kd1 ka kd qm2 k12 k21 ka2 kd2 beta s v global NI QM KA1 KD1 KA KD QM2 K12 K21 KA2 KD2 BETA S V if strcmp(it,'sma') V = repmat(v(2:end) , S = repmat(s(2:end) , KA1 = repmat(ka1(2:end), KD1 = repmat(kd1(2:end), end  [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]);  if or(strcmp(it,'scl'),strcmp(it,'mcl')) NI = repmat(repmat(ni,1,n), [zs,1,nc]); QM = repmat(qm , [zs,1,nc]); KA = repmat(ka , [zs,1,nc]); KD = repmat(kd , [zs,1,nc]); end  166  if strcmp(it,'bil') NI = repmat(repmat(ni,1,n), QM = repmat(qm , KA = repmat(ka , KD = repmat(kd , QM2 = repmat(qm2 , KA2 = repmat(ka2 , KD2 = repmat(kd2 , end  [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]);  if strcmp(it,'spr') NI = repmat(repmat(ni,1,n), QM = repmat(qm , KA = repmat(ka , KD = repmat(kd , BETA = repmat(beta , K12 = repmat(k12 , K21 = repmat(k21 , KA2 = repmat(ka2 , KD2 = repmat(kd2 , end  [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]);  if strcmp(it,'sprS') NI = repmat(repmat(ni,1,n), QM = repmat(qm , KA = repmat(ka , KD = repmat(kd , BETA = repmat(beta , K12 = repmat(k12 , K21 = repmat(k21 , End  [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]); [zs,1,nc]);  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function blockdiag is called by “discrete”. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function C = blockdiag(A,B) a = size(A); b = size(B); C = [A spalloc(a(1),b(2),0) ; spalloc(b(1),a(2),0) B];  167  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “residual_zonal” is called by “optimize_zonal”; it % simulates a breakthrough curve using the inputted parameters and then % calculates, and returns to “optimize_zonal”, the residual between % this calculated BTC and a given BTC. % In this version of the function the parameters that are being % optimized are the isotherm parameters ‘ka’ ‘kd’ and ‘qm’ for four % flowrates: 1.5, 5, 10 and 20 mL/min. ‘ka’ and ‘kd’ are flowrate % independent whereas ‘qm’ varies with flow. Only experimental data % under binding conditions is being used for parameter optimization. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function res = residual_zonal(param_log) global t t0 te nc nt nd ec sv global zs ts global data  % necessary % derivable % defined  global zr qstar qm qm2 beta ka ka1 ka2 kd kd1 kd2 k12 k21 Kf d s v Ez u col_area param = 10.^param_log;  % translates the parameters from % “optimize_zonal” out of logarithmic % scale  res = []; for i = 1:nd if data(i).comp switch i case {1}; uu = 1.5; zr = [125.32 61.69 171.57]; timelag = 15.245; qm = param(3); case {2}; uu = 5; zr = [37.60 27.93 77.67]; timelag = 4.5735; qm = param(4); case {3}; uu = 10; zr = [37.60 27.93 77.67]; timelag = 2.2868; qm = param(5); case {4}; uu = 20; zr = [9.40 8.35 23.22] timelag = 1.1434; qm = param(6); end u  % uu is flow in mL/min  = uu * (1/60) * (1/100^3) * (1/col_area) * (1/ec);  168  ka = param(1); kd = param(2); %ensures that there is a timepoint corresponding to the timelag te = data(i).tm(end); if not(timelag < te); data(i).res = Inf*data(i).Cm; break; end if (timelag) > 0; t0 = [timelag te]; else t0 = te; end discrete;  % calls “discrete” function with new parameter values  % time vector t  = data(i).tm;  j = i + 1; while not(j > nd) if data(j).comp == 0; t = union(t,data(j).tm); else j = nd; end j = j + 1; end ts = length(t); if t0(1) < t(2); t = union(t0(1)/2,t) ; end solve_zonal;  % solving occurs here  if zm == 10 Cs = C(:,zs,end); else Cs = Ct(:,nt,end); end end % residual data(i).res = data(i).Cm - Cs(ismember(t,data(i).tm)); switch i case {1,2,3,4}; if out; figure(1); hold on; plot(data(i).tm,data(i).Cm,'k.'); plot(t,Cs,'r'); drawnow; end end res = [res ; data(i).res]; end  169  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “solve_zonal” is the main solution program; it returns a % BTC using inputted parameter values. The functions “initial_zonal”, % “state_zonal_to” and “state_zonal_from” are called in “solve_zonal”. % The built-in function ode15s is used to solve the equations defined % in “right_zonal”. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function solve_zonal global t at rt dbl dbr t0 te n it qs zi nc nt % necessary global zs ts % derivable global Ct C Ci Q Q2 % defined % Initial conditions & mass matrix if strcmp(it,'scl') || strcmp(it,'mcl') || strcmp(it,'sma') [ct0 c0 ci0 q0] = initial_zonal; elseif strcmp(it,'bil') || strcmp(it,'spr') || strcmp(it,'sprS') [ct0 c0 ci0 q0 q20] = initial_zonal2; end mmt = ones(n,nt); if strcmp(it,'scl') || strcmp(it,'mcl') || strcmp(it,'sma') mmc = ones(zs,n,nc,3); if strcmp(it,'sma') mmc(:,1,:,3) = 0; end if qs mmc(:,:,:,3) = 0; end elseif strcmp(it,'bil') || strcmp(it,'spr') || strcmp(it,'sprS') mmc = ones(zs,n,nc,4); if qs mmc(:,:,:,3) = 0; mmc(:,:,:,4) = 0; end end if dbl mmc(1,:,:,1) = 0; end if dbr mmc(zs,:,:,1) = 0; end mm = [mmt(:) ; mmc(:)]; if strcmp(it,'scl') || strcmp(it,'mcl') || strcmp(it,'sma') y0 = state_zonal_to(ct0,c0,ci0,q0);  170  elseif strcmp(it,'bil') || strcmp(it,'spr') || strcmp(it,'sprS') y0 = state_zonal_to2(ct0,c0,ci0,q0,q20); end MM = diag(sparse(mm)); if t0(end) < te zn = length(t0) + 1; else zn = length(t0); end if strcmp(it,'scl') || strcmp(it,'mcl') || strcmp(it,'sma') if zn > 1 for i = 1:zn-1 if isempty(find(t == t0(i),1)); t = sort([t t0(i)]); end; p(i) = find(t==t0(i),1); end p = [p length(t)]; [T Y] = ode15s(@right_zonal,t(1:p(1)),y0,… odeset('Stats','off','Mass',MM,'AbsTol',at,'RelTol',rt),1); for i = 2:zn y0 = Y(end,:)'; [th Yh] = ode15s(@right_zonal,t(p(i-1):p(i)),y0,… odeset('Stats','off','Mass',MM,'AbsTol',at,'RelTol',… rt),i); Y = [Y ; Yh(2:end,:)]; end else [t Y] = ode15s(@right_zonal,t,y0,odeset('Stats','off','Mass',… MM,'AbsTol',at,'RelTol',rt),1); end % Reshaping of solution matrix [Ct C Ci Q] = state_zonal_from(Y); elseif strcmp(it,'bil') || strcmp(it,'spr') || strcmp(it,'sprS') if zn > 1 for i = 1:zn-1 if isempty(find(t == t0(i),1)); t = sort([t t0(i)]); end; p(i) = find(t == t0(i),1); end p = [p length(t)]; [T Y] = ode15s(@right_zonal,t(1:p(1)),y0,… odeset('Stats','off','Mass',MM,'AbsTol',at,'RelTol',rt),1); for i=2:zn y0 = Y(end,:)'; [th Yh] = ode15s(@right_zonal,t(p(i-1):p(i)),y0,… odeset('Stats','off','Mass',MM,'AbsTol',at,'RelTol',… rt),i); Y = [Y ; Yh(2:end,:)]; end  171  else [t Y] = ode15s(@right_zonal,t,y0,… odeset('Stats','off','Mass',MM,'AbsTol',at,'RelTol',… rt),1); end % Reshaping of solution matrix [Ct C Ci Q Q2] = state_zonal_from2(Y); end % Inlet profile if not(dbl) for i = 1:ts for j = 1:nc C(i,1,:,j) = Ct(i,zi(j),:); end end end  172  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The functions “initial_zonal” and “initial_zonal2” define the initial % concentrations in the column. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [ct c ci q] = initial_zonal global zs n it nc nt ct c ci q  = = = =  zeros(nt,n); zeros(zs,n,nc); zeros(zs,n,nc); zeros(zs,n,nc);  % Steric Mass Action (SMA) global c0 d if strcmp(it,'sma') ct(:,1) = c0(1); c(:,1,:) = c0(1); ci(:,1,:) = c0(1); q(:,1,:) = d; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [ct c ci q q2] = initial_zonal2 global zs n nc nt ct c ci q q2  = = = = =  zeros(nt,n); zeros(zs,n,nc); zeros(zs,n,nc); zeros(zs,n,nc); zeros(zs,n,nc);  173  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “right_zonal” defines the equations to be solved for % every configuration of the Zonal Rate Model as well as the boundary % conditions. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dy = right_zonal(t,y,zn) global zs M2 it EZ KF Ez u dbl dbr fd n ep zm nc zi zr % Memory allocation c1 = zeros(zs,n,nc); c2 = zeros(zs,n,nc); % Reshaping of state vector if strcmp(it,'scl') || strcmp(it,'mcl') || strcmp(it,'sma') [ct c ci q] = state_zonal_from(y); elseif strcmp(it,'bil') || strcmp(it,'spr') || strcmp(it,'sprS') [ct c ci q q2] = state_zonal_from2(y); end % CSTR network equations switch zm case 10; dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); case 11; dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(1) * (c(zs,:,1) - ct(2,:)); case 13; dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (c(zs,:,1) - ct(2,:)); case 14; dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(2) * (c(zs,:,1) - ct(3,:)); dct(4,:) = 1/zr(3) * (ct(3,:) - ct(4,:)); case 15; dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(3) * (c(zs,:,1) - ct(3,:)); dct(4,:) = 1/zr(4) * (ct(3,:) - ct(4,:)); case 16; dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(3) * (ct(2,:) - ct(3,:)); dct(4,:) = 1/zr(4) * (c(zs,:,1) - ct(4,:)); dct(5,:) = 1/zr(5) * (ct(4,:) - ct(5,:)); dct(6,:) = 1/zr(6) * (ct(5,:) - ct(6,:)); case 21; dct(1,:) = (1/zr(2)+1/zr(3)) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(1) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(1) * (c(zs,:,2) - ct(3,:)); dct(4,:) = 1/zr(2) * (c(zs,:,1) - ct(4,:)) + 1/zr(3) *… (ct(3,:) - ct(4,:)); case 22; dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(2) * (c(zs,:,2) - ct(3,:));  174  case 23;  case 24;  case 25;  case 26;  case 31;  case 33;  case 34;  case 36;  dct(4,:) = 1/zr(3) * (c(zs,:,1) - ct(4,:)) + 1/zr(4) *… (ct(3,:) - ct(4,:)); dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(3) * (c(zs,:,2) - ct(3,:)); dct(4,:) = 1/zr(4) * (c(zs,:,1) - ct(4,:)) + 1/zr(5) *… (ct(3,:) - ct(4,:)); dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(3) * (ct(1,:) - ct(3,:)); dct(4,:) = 1/zr(2) * (c(zs,:,1) - ct(4,:)); dct(5,:) = 1/zr(3) * (c(zs,:,2) - ct(5,:)); dct(6,:) = 1/zr(4) * (ct(4,:) - ct(6,:)) + 1/zr(5) *… (ct(5,:) - ct(6,:)); dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(3) * (ct(1,:) - ct(3,:)); dct(4,:) = 1/zr(4) * (c(zs,:,1) - ct(4,:)); dct(5,:) = 1/zr(5) * (c(zs,:,2) - ct(5,:)); dct(6,:) = 1/zr(6) * (ct(4,:) - ct(6,:)) + 1/zr(7) *… (ct(5,:) - ct(6,:)); dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(3) * (c(zs,:,1) - ct(3,:)); dct(4,:) = 1/zr(4) * (c(zs,:,2) - ct(4,:)) + 1/zr(5) *… (ct(3,:) - ct(4,:)); dct(1,:) = (1/zr(4)+1/zr(5)) * (inlet(t,zn) - ct(1,:)); dct(2,:) = (1/zr(2)+1/zr(3)) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(1) * (ct(2,:) - ct(3,:)); dct(4,:) = 1/zr(1) * (c(zs,:,3) - ct(4,:)); dct(5,:) = 1/zr(2) * (c(zs,:,2) - ct(5,:)) + 1/zr(3) *… (ct(4,:) - ct(5,:)); dct(6,:) = 1/zr(4) * (c(zs,:,1) - ct(6,:)) + 1/zr(5) *… (ct(5,:) - ct(6,:)); dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(3) * (ct(2,:) - ct(3,:)); dct(4,:) = 1/zr(4) * (c(zs,:,3) - ct(4,:)); dct(5,:) = 1/zr(5) * (c(zs,:,2) - ct(5,:)) + 1/zr(6) *… (ct(4,:) - ct(5,:)); dct(6,:) = 1/zr(7) * (c(zs,:,1) - ct(6,:)) + 1/zr(8) *… (ct(5,:) - ct(6,:)); dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(3) * (ct(1,:) - ct(3,:)); dct(4,:) = 1/zr(4) * (ct(1,:) - ct(4,:)); dct(5,:) = 1/zr(2) * (c(zs,:,1) - ct(5,:)); dct(6,:) = 1/zr(3) * (c(zs,:,2) - ct(6,:)); dct(7,:) = 1/zr(4) * (c(zs,:,3) - ct(7,:)); dct(8,:) = 1/zr(5) * (ct(5,:) - ct(8,:)) + 1/zr(6) *… (ct(6,:) - ct(8,:)) + 1/zr(7) * (ct(7,:) - ct(8,:)); dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(3) * (ct(2,:) - ct(3,:)); dct(4,:) = 1/zr(4) * (c(zs,:,1) - ct(4,:)); dct(5,:) = 1/zr(5) * (c(zs,:,2) - ct(5,:)) + 1/zr(6) *… (ct(4,:) - ct(5,:));  175  dct(6,:) = 1/zr(7) * (c(zs,:,3) - ct(6,:)) + 1/zr(8) *… (ct(5,:) - ct(6,:)); case 44; dct(1,:) = 1/zr(1) * (inlet(t,zn) - ct(1,:)); dct(2,:) = 1/zr(2) * (ct(1,:) - ct(2,:)); dct(3,:) = 1/zr(3) * (ct(1,:) - ct(3,:)); dct(4,:) = 1/zr(4) * (ct(1,:) - ct(4,:)); dct(5,:) = 1/zr(5) * (ct(1,:) - ct(5,:)); dct(6,:) = 1/zr(2) * (c(zs,:,1) - ct(6,:)); dct(7,:) = 1/zr(3) * (c(zs,:,2) - ct(7,:)); dct(8,:) = 1/zr(4) * (c(zs,:,3) - ct(8,:)); dct(9,:) = 1/zr(5) * (c(zs,:,4) - ct(9,:)); dct(10,:) = 1/zr(6) * (ct(6,:) - ct(10,:)) + 1/zr(7) *… (ct(7,:) - ct(10,:)) + 1/zr(8) * (ct(8,:) - ct(10,:)) + … 1/zr(9) * (ct(9,:) - ct(10,:)); end % Left membrane boundary if not(dbl) for i = 1:nc c(1,:,i) = ct(zi(i),:); end end for i = 1:nc c1(:,:,i) = first_derivative(c(:,:,i)); c2(:,:,i) = M2*c(:,:,i); end % Continuity equation if fd == 0 dg = -u*c1 + EZ.*c2; dh = zeros(zs,n,nc); if strcmp(it,'scl') || strcmp(it,'mcl') || strcmp(it,'sma') dq = isotherm(c,q); elseif strcmp(it,'bil') [dq dq2] = isothermbil(c,q,q2); elseif strcmp(it,'spr') [dq dq2] = isothermspr(c,q,q2); elseif strcmp(it,'sprS') [dq dq2] = isothermsprS(c,q,q2); end else dg = -u*c1 + EZ.*c2; dh = 1/ep*KF.*(c - ci); if strcmp(it,'scl') || strcmp(it,'mcl') || strcmp(it,'sma') dq = isotherm(ci,q); elseif strcmp(it,'bil') [dq dq2] = isothermbil(ci,q,q2); elseif strcmp(it,'spr') [dq dq2] = isothermspr(c,q,q2); elseif strcmp(it,'sprS') [dq dq2] = isothermsprS(c,q,q2); end end  176  % Membrane boundaries if dbl for i = 1:nc dg(1,:,i) = Ez.*c1(1,:,i) - u*(c(1,:,i) - ct(zi(i),:)); end else dg(1,:,:) = 0; end if dbr dg(zs,:,:) = c1(zs,:,:); end % Right hand side of differential-algebraic equation system if strcmp(it,'scl') || strcmp(it,'mcl') || strcmp(it,'sma') dy = [dct(:) ; dg(:) ; dh(:) ; dq(:)]; elseif strcmp(it,'bil') || strcmp(it,'spr') || strcmp(it,'sprS') dy = [dct(:) ; dg(:) ; dh(:) ; dq(:) ; dq2(:)]; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “firstderivative” is called above. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function c1 = first_derivative(c) global A0 c1 = A0*c; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “inlet” sets the inlet concentration. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function c = inlet(t,zn) global c0 t0 it if length(t0) > 1 c = c0 * (zn == 2); else c = c0 * (zn == 1); end if strcmp(it,'sma') c(1) = c0(1); end  177  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “isotherm” defines the Langmuir and SMA isotherms. % The function “isothermbil” defines the bi-Langmuir isotherm. % The function “isothermspr” defines the full spreading isotherm. % The function “isothermsprS” defines the simplified spreading % isotherm. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dq = isotherm(c,q) global it % Steric Mass Action (SMA) global V S KA1 KD1 d n if strcmp(it,'sma') cp = c(:,2:end,:); qp = q(:,2:end,:); qb = q(:,1,:) - sum(S.*qp,2); % Ions available for binding cs = repmat(c(:,1,:),1,n-1); qb = repmat(qb,1,n-1); dq(:,1,:) = d - q(:,1,:) - sum(V.*qp,2); % Electroneutrality dq(:,2:n,:) = KA1.*cp.*qb.^V - KD1.*qp.*cs.^V; % SMA Kinetic end % Langmuir global KA KD QM NI if strcmp(it,'scl') dq = KA.*c.*(QM - NI.*q) - KD.*q; end if strcmp(it,'mcl') dq = KA.*c.*QM.*(1 - NI.*repmat(sum(q./QM,2),1,n)) - KD.*q; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [dq dq2] = isothermbil(c,q,q2) global KA KD QM NI KA2 KD2 QM2 dq = KA.*c.*(QM - NI.*q) - KD.*q; dq2 = KA2.*c.*(QM2 - NI.*q2) - KD2.*q2;  178  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [dq dq2] = isothermspr(c,q,q2) global KA KD QM NI K12 K21 KA2 KD2 BETA dq = (KA.*c - K12.*q).*(QM - NI.*q - NI.*BETA.*q2) - KD.*q + K21.*q2; dq2 = (KA2.*c + K12.*q).*(QM - NI.*q - NI.*BETA.*q2) - KD2.*q2 -… K21.*q2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [dq dq2] = isothermsprS(c,q,q2) global KA KD QM NI K12 K21 BETA dq = (KA.*c - K12.*q).*(QM - NI.*q - NI.*BETA.*q2) - KD.*q + K21.*q2; dq2 = K12.*q.*(QM - NI.*q - NI.*BETA.*q2) - K21.*q2;  179  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The functions “state_zonal_to” and “state_zonal_to2” define the % continuity equation for membrane chromatography. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function y = state_zonal_to(ct,c,ci,q) global fd ec ep zs n nc switch fd case 0 g = c + (1 - ec)/ec * q; h = zeros(zs,n,nc); case 1 g = c + (1 - ec)/ec * (ep*ci + (1 - ep)*q); h = ci + (1 - ep)/ep * q; end y = [ct(:) ; g(:) ; h(:) ; q(:)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function y = state_zonal_to2(ct,c,ci,q,q2) global fd ec ep zs n nc switch fd case 0 g = c + (1 - ec)/ec * (q + q2); h = zeros(zs,n,nc); case 1 g = c + (1 - ec)/ec * (ep*ci + (1 - ep)*(q + q2)); h = ci + (1 - ep)/ep * (q + q2); end y = [ct(:) ; g(:) ; h(:) ; q(:) ; q2(:)];  180  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The functions “state_zonal_from” and “state_zonal_from2” reshape the % solution matrix. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [ct c ci q] = state_zonal_from(y) global fd ec ep zs n nt nc if size(y,2) == 1 ct = reshape(y(1:nt*n),nt,n); y = reshape(y(nt*n+1:end),zs,n,nc,3); g = y(:,:,:,1); h = y(:,:,:,2); q = y(:,:,:,3); else ct = reshape(y(:,1:nt*n),[],nt,n); y = reshape(y(:,n*nt+1:end),[],zs,n,nc,3); g = y(:,:,:,:,1); h = y(:,:,:,:,2); q = y(:,:,:,:,3); end switch case 0 c = ci = case 1 c = ci = end  fd g - (1 - ec)/ec * q; zeros(zs,n); g - (1 - ec)/ec * ep * h; h - (1 - ep)/ep * q;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [ct c ci q q2] = state_zonal_from2(y) global fd ec ep zs n nt nc if size(y,2) == 1 ct = reshape(y(1:nt*n),nt,n); y = reshape(y(nt*n+1:end),zs,n,nc,4); g = y(:,:,:,1); h = y(:,:,:,2); q = y(:,:,:,3); q2 = y(:,:,:,4); else ct = reshape(y(:,1:nt*n),[],nt,n); y = reshape(y(:,n*nt+1:end),[],zs,n,nc,4); g = y(:,:,:,:,1); h = y(:,:,:,:,2); q = y(:,:,:,:,3); q2 = y(:,:,:,:,4);  181  end switch case 0 c = ci = case 1 c = ci = end  fd g - (1 - ec)/ec * (q + q2); zeros(zs,n); g - (1 - ec)/ec * ep * h; h - (1 - ep)/ep * (q + q2);  182  Appendix F: Program code for the moment analysis of breakthrough curves The first absolute and second central moment of a BTC are given by [207]:   1   dcout tdt  0 0 dt  2   dcout t  1 2 dt  0 0 dt  1  1    (F-1)      (F-2)  where, 0, the zeroth moment is defined as:   0   0  dcout dt dt  (F-3)  The following MATLAB program calculates a BTC using the ZRM then determines its first absolute and second central moments using the above equations and a numerically derived derivative of the BTC. The integrals are calculated using the built-in numerical integrator trapz which uses the trapezoidal method of integration.  183  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % This program uses the ZRM to calculate a BTC for a given set of % parameters then uses the built-in numerical integrator trapz to % calculate moments based on a numerical differentiation of the BTC. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function mom_surf clear global global t t0 te timelag ts global Ct Cs tt zr col_area Ez u ze global first_moment second_moment BT_C time parameter; Pe = [100 1000 10000 50000]; % range of Pe number F = [0:0.1:1]; % flow fraction through inner zone i.e. Q1/Q ts = 400; % time steps r1 = 151.75; r4 = 49.41; timelag = 17.17;  % model parameters  tend = [2500 2500 2500 2500]; %u is interstitial velocity (m/s) u = 1.5 * (1/60) * (1/100^3) * (1/col_area) * (1/ec); ze = 2.2e-3;  % column length m  for i = 1:length(Pe) Ez = ze*u/Pe(i); te = tend(i); t0 = [timelag te];  % axial dispersion coeff. % calculation end time  for j = 1:length(F) zr = [r1 r4/(1 - F(j)) r4/F(j)]; % with given zr and parameter values the ZRM is used to produce % a breakthrough curve whose moments can then be analyzed discrete; if t0(1)<t(2); t = union(t0(1)/2,t) ; end solve_zonal;  % ZRM produces BTC  tt = t';  184  Cs = Ct(:,nt,1); % numerically determine dC/dt for k = 1:length(Cs)-1 dCsdt(k,1) = (Cs(k+1)-Cs(k))/(tt(k+1)-tt(k)); end % built in integrator determines moments mom1 = trapz(tt(1:end-1),dCsdt.*tt(1:end-1)); mom0 = trapz(tt(1:end-1),dCsdt); mom2 = trapz(tt(1:end-1),dCsdt.*(tt(1:end-1)-mom1/mom0).^2); first_moment(i,j) = mom1/mom0; second_moment(i,j) = mom2/mom0;  % first absolute moment % second central moment  BT_C(:,(j+(i-1)*length(F))) = Cs; % concentration for BTC time(:,(j+(i-1)*length(F))) = tt; % time for BTC end end  185  Appendix G: Program code for the determination of Thomas model parameters The following MATLAB programs optimize the Langmuir isotherm parameters ka and kd by fitting experimental BTC data to results from the model of membrane chromatography described in Section 5.4.2 that consists of the analytical Thomas model sandwiched between two CSTRs. To produce a model-derived BTC the program calculates the solute concentration profile at the outlet to the first CSTR using Equation 4-1. It then uses this profile as the inlet to the Thomas model (Equations 5-18 through 5-22) whose resultant ouptut is then used as input to a second CSTR.  The output of the second CSTR  represents the model-derived BTC.  186  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “optimize_Thomas” similarly to “optimize_zonal”, % optimizes the Langmuir isotherm parameters ‘ka’ and ‘kd’ in order to % fit a given BTC using the Thomas model sandwiched between two CSTRs. % It minimizes the residual calculated by “residual_Thomas” using the % built-in function lsqonlin and then returns the optimized parameter % values. % This version of the program optimizes parameters for binding % conditions at 20 mL/min. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function optimize_Thomas clear global; global data cin e L qm ka kd nd col_area uu tt % Parameters e = 0.7; L = 2.2e-3; qm = 402; cin = 1; col_area = 2.206e-3; uu = 20; tt = [0:10:2500]; data(1).file data(1).comp  % flowrate in mL/min % time domain  = 'avg-binding-20mLmin-1mgmL.mat'; = 1;  % Loads data load(data(1).file); data(1).tm = tm; data(1).Cm = Cm; param = [1.12e-2 1.11e-3]; param_log = log10(param);  % Logarithmic scale  sum(residual_Thomas(param_log).^2) param_opt_log = lsqnonlin(@residual_Thomas,param_log,… 0*param_log-16,0*param_log+16,optimset('display','iter',… 'DiffMinChange',1e-4,'MaxFunEvals',Inf,'TolX',1e-4)); dev_opt = residual_Thomas(param_opt_log); res_lsq = sum(dev_opt.^2); param_opt = 10.^param_opt_log param_opt_log res_lsq  187  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “residual_Thomas” is called by “optimize_Thomas”; it % simulates a breakthrough curve using the inputted parameters and then % calculates, and returns to “optimize_Thomas”, the residual between % this calculated BTC and a given BTC. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function res = residual_Thomas(param_log) global data nd cin e L qm ka kd col_area v uu tlag tau tt t c param = 10.^param_log;  % translates out of Logarithmic scale  figure(1); clf res = []; for i = 1:nd ka = param(1); kd = param(2); v = uu * (1/60) * (1/100^3) * (1/col_area) * (1/e); tlag = 0.38/uu*60; % determined from ZRM work tau = 1/(4.5716e-3 + 0.0053565*uu); % under flowthrough conditions Thomas_model_func;  % solving occurs here  dumt = data(i).tm; dumcm = data(i).Cm; dumc = interp1(t,c,dumt)';  data(i).res = dumcm - dumc;  % interpolates calculated BTC to % get concentrations at the same % time points as the data set % calculates the residual  switch i case {1,2,3,4}; figure(1); hold on; plot(data(i).tm,data(i).Cm,'k.'); plot(t,c,'r'); drawnow; end res = [res ; data(i).res]; end  188  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “Thomas_model_func” produces a BTC based on the % analytical solution of the Thomas model sandwiched between two CSTRs. % This program solves the Thomas model using the built in functions % besseli which yields the zeroth order Bessel Function and trapz which % numerically evaluates an integer. % The differential equation representing the downstream CSTR is solved % using the built-in function ode15s which calls the function % “Thomas_tank2” where the equation is defined. % This version of the program optimizes parameters for binding % conditions at 20 mL/min. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Thomas_model_func global tau e L qm cin v ka kd tt z CC Kd n t c uu tlag z = length(tt); Kd = kd/ka; n = (1 - e)*qm*ka*L/e/v; r = zeros(1,z); T = zeros(1,z); J1 = zeros(1,z); J2 = zeros(1,z); for i = 1:z r(i) = 1 + (cin - cin*exp(-tt(i)/tau))/Kd; % solves the first CSTR T(i) = e*Kd*r(i)/(1 - e)/qm*(v*tt(i)/L - 1); int = @(x) exp(-x).*besseli(0,2*(x.*(n*T(i))).^0.5); J1(i) = 1 - exp(-n*T(i))*trapz(int,0,n/r(i)); int2 = @(x) exp(-x).*besseli(0,2*(x.*(n*T(i)/r(i))).^0.5); J2(i) = 1 - exp(-n*T(i)/r(i))*trapz(int2,0,n); C(i) = J1(i)/(J1(i) + (1 - J2(i))*exp((1 - 1/r(i))*(n -n*T(i)))); end CC = real(C); c0 = 0; % must ensure that the inlet to the second tank occurs at an % appropriately zeroed time, hence the 60 seconds. tt = tt - 60; tt = tt(7:end); CC = CC(7:end); eqn = @(t,y) Thomas_tank2(t,y,tau,tt,CC); [t c] = ode15s(eqn,[0:5:2500],c0); t = t + 60 + tlag;  189  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % The function “Thomas_tank2” defines the differential equation for the % downstream CSTR. It uses the built in function interp1 to % interpolate values from the tank inlet profile passed from % “Thomas_model_func”. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dy = Thomas_tank2(t, y, tau, tt, CC) dy = 1/tau*(interp1(tt,CC,t) - y);  190  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0059076/manifest

Comment

Related Items