PROTEIN POLARIZATION IN PACKED HOLLOW FIBRE BIOREACTORS by JURGEN KOSKA Dipl. Ing., Mannheim Institute of Technology, Germany, 1991 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUEST FOR THE DEGREE OF MASTER IN APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Chemical Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA August 1993 ©JurgenKoska,193 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of (%-eee( ice," Al; ‘ z The University of British Columbia Vancouver, Canada Date DE-6 (2/88) ^r, - Abstract Osmotically active proteins and cells, retained in the extracapillary space (ECS) of hollow fibre bioreactors (HFBRs), can influence the hydrodynamics of such devices. A mathematical model was developed to describe the coupled hydrodynamics and high molecular weight protein transport in a cell filled HFBR. It was assumed that the multi-fibre reactor can be represented by a single, straight fibre surrounded by a symmetry envelope containing fluid and a homogeneous packed bed of cells. The low Reynolds number flow in this porous medium was described by Darcy's law. Because of difficulties associated with operating a reactor filled with mammalian cells, a suitable analogue was used to experimentally investigate protein polarization in packed HFBRs. The ECS side of the reactor was filled with an agarose/protein solution which, upon cooling, formed a porous medium with a uniform initial concentration of protein. A constant lumen flow was established for several days before the cartridge was sacrificed and the axial ECS protein distribution was measured. Since the protein transport and HFBR hydrodynamics were coupled, numerical methods were required to solve the governing equations of both the two-dimensional (axial plus radial variations) and the one-dimensional (axial variations only) models developed to predict axisymmetric transient ECS protein concentrations. Computer modelling results indicated that, because of the large length/radius ratio of the representative fibre unit, the twodimensional ECS protein concentrations could be accurately duplicated by the simpler onedimensional model. The latter model, required about two orders-of-magnitude less computational time than the former. The one-dimensional model results were compared to experimentally obtained ECS protein profiles and, subsequently, the model was used to predict protein polarization in ITFBRs for different conditions. The hydraulic conductivity of the agarose gel, required for the model, was experimentally determined using the falling head method. The measured conductivity values failed to iii adequately describe the observed protein polarization in the ECS of HFBRs. However, by using a gel conductivity which was about an order-of-magnitude higher than the measured value, it was found that the model agreed well, in general, with the measured ECS protein polarization profiles obtained for initial protein loadings of 5 - 30 g/L. The higher apparent conductivity, needed to fit the model to the experimental results, was attributed to the inability of the gel to completely fill the geometrically complex ECS. Since in HFBRs, packed anchorage-dependent mammalian cells are expected to achieve hydraulic conductivities similar to values encountered in tissues, several model simulations were carried out at low tissue conductivity values. The results indicated that, for these conditions, protein transport in the ECS is mainly governed by diffusion. Protein polarization, a dominant feature of empty ECS protein transport, is greatly reduced. Also, it was shown that the removal of product protein from a packed ECS space can be difficult, since ECS flows are reduced. Models, such as those developed here, can be used to further investigate HFBR operation and process optimization. iv Table of Contents Abstract ^ Table of Contents ^ iv List of Tables ^ vi List of Figures ^ vii Acknowledgments ^ CHAPTER 1 Introduction ^ CHAPTER 2 Literature Review ^ 2.1. Cell Culture ^ 2.2. Hollow Fibre Bioreactor Modelling ^ 2.3. Protein Transport in a Gel ^ 2.4. Porous Bed Conductivity ^ CHAPTER 3 Model Formulation ^ 3.1. Krogh Cylinder Approximation ^ 3.2. Hydrodynamics ^ 3.2.1. Flow in Porous Media ^ 3.2.2. Membrane Hydrodynamics ^ 3.2.3. Lumen Hydrodynamics ^ 3.2.4. Extracapillary Space Hydrodynamics ^ 3.3. Protein Transport ^ 3.3.1. Protein Properties ^ 3.3.2. Governing Equations ^ 3.3.3. Simplified One-Dimensional Model ^ 3.4. Numerical Methods ^ CHAPTER 4 Materials and Methods ^ 4.1. Hollow Fibre Bioreactor Specifications ^ 4.1.1. Fibre Dimensions ^ 4.1.2. Membrane Permeability Measurements ^ 1 5 5 6 10 11 14 14 15 15 17 19 22 23 23 27 28 30 33 33 33 34 4.2. Gel Preparation ^ 4.3. Reactor Preparation and Flow Configuration ^ 4.3.1. Preparation and Analysis of Samples ^ 4.4. Gel Conductivity Measurements ^ 35 36 37 39 CHAPTER 5 Results and Discussion ^ 41 5.1. Reactor Specifications ^ 41 5.2. Gel Conductivity Measurements ^ 45 5.3. Model Parameters ^ 47 5.4. Modelling of Hollow Fibre Bioreactor Hydrodynamics ^ 50 5.5. Protein Transport in a Packed Hollow Fibre Bioreactor ^ 54 5.5.1. Comparison of Experimental Data with Model Results ^ 56 62 5.5.2. Model Predictions ^ CHAPTER 6 Conclusions ^ 6.1. Experimental Versus Model Results ^ 76 76 vi List of Tables 2.1 Conductivities of packed spheres ^ 12 2.2 Hydraulic bed conductivities ^ 13 3.1 Influence of BSA concentration on viscosity, density and diffusivity; percent change in value from an infinitely dilute solution ^ 26 5.1 Summary of gambro ® GFE 15 reactor specifications ^ 44 5.2 Experimentally determined agarose gel conductivities ^ 46 5.3 Modelling parameters for HFBR protein polarization program ^ 49 5.4 Summary of experimental protein polarization results ^ 57 5.5 Maximum average axial ECS Peclet numbers ^ 73 vii List of Figures 1.1^Pressure distribution along a closed-shell HFBR ^ 2 1.2^Starling flow around a single fibre in a closed-shell HFBR ^ 2 3.1^Single fibre Krogh cylinder geometry ^ 14 3.2^Schematic of the Krogh cylinder representation of a multi-fibre reactor ^ 15 3.3^Osmotic pressure of BSA protein solutions. Dashed line: Vilker's correlation (Equation 3.23),.: Scatchard's data, solid line: fit of Equation. 3.22 to Scatchard's data ^ 25 3.4 Computational control volume ^ 32 4.1^Schematic diagram of the reactor flow configuration ^ 37 4.2 Schematic diagram of the falling head method for gel conductivity measurements ^ 39 5.1 Measured dimensions of the gambro® GFE 15 HFBR cartridge ^ 43 5.2^Axial ECS pressures at R2 calculated using Apelblat's analytical solution (-) and the 2-D (.) and 1-D (0) computer models. k3= 5x10-15 m2, 6 = 1.0, API = 5000 Pa, Lp= 6x10-15m ^ 50 5.3 Axial ECS pressure distribution for different conductivities k3, Q 500 mL/min, c = 1.0, Lp= 6x10-15 m ^ 51 5.4 Transmembrane and average ECS velocities for different conductivities k3. Q = 500 mL/min, = 1.0, Lp= 6x10'5 m. Solid line: transmembrane velocity, dotted line: radially-averaged axial ECS velocity 52 5.5^Axial ECS pressure distribution for different membrane hydraulic permeabilities k3 = 5x10-15 m2, 6 = 1.0, Q = 500 mL/min ^ 53 viii 5.6 Transient axial protein concentration results from 2-D (-) and 1-D (.) model. The 2-D model results represent concentrations at R2 and R3, whereas the 1-D results represent radially-averaged concentrations c30 = 10 g/L, Q= 500 mL/min, k3 = 5x10-15 M2, = 1.0, Lp= 6X10-15m, 55 5.7 Experimentally obtained radially-averaged axial ECS protein concentrations versus computer modelling results. c30= 20.2 g/L, 0 = 492 mL/min, time = 215.6 h, 3% gel, k3= 5x10-15 and 7.5x10-16 m2, 6 = 1.0, Lp= 6x10-15m. Solid line: model; •: experimental data, dashed line: section length 59 5.8 Experimentally obtained radially-averaged axial ECS protein concentrations versus computer modelling results. c30 = 19.3 g/L, Q= 544 mL/min, time = 116.3 h, 3% gel, k3= 5x10-15 m2, c = 1.0, Lp= 6x10-15m. Solid line: model;.: experimental data 59 5.9 Experimentally obtained radially-averaged axial ECS protein concentrations versus computer modelling results. c30= 27.5 g/L, Q= 512 mL/min, time = 331.5 h, 2% gel, k3= 5X10-15M2, c = 1.0, Lp= 6X10-15m. Solid line: model;.: experimental data 61 5.10 Experimentally obtained radially-averaged axial ECS protein concentrations versus computer modelling results. c30= 26.9 g/L, Q= 447 mL/min, time = 380.8 h, 2% gel, k3= 5X10-15M2, c = 1.0, Lp= 6X10-15 m. Solid line: model; dotted line: steadystate result at 1888 h; •: experimental data 61 5.11 Experimentally obtained radially-averaged axial ECS protein concentrations versus computer modelling results. c30 = 4.7 g/L, Q= 484 mL/min, time = 329.2 h, 2% gel, k3= 5x10-15 m2, 6 = 1.0, Lp= 6X10-15m. Solid line: model;.: experimental data 62 5.12 Lumen and ECS total pressures (P3 - 113) as function of axial distance at steadystate for a packed ECS conductivity. c30= 10 g/L, Q= 500 mL/min, k3 = 10-13 m2, 6 = 1.0,Lp= 6x10-15m 63 5.13 Lumen and ECS total pressures (P3- H3) as function of axial distance at steadystate for an empty ECS. c30= 10 g/L, Q= 500 mL/min, k3 = 10-9 m2, Lp= 6x10-15m 63 5.14 Transient radially-averaged axial ECS protein distribution and axial velocities for an empty ECS. c30 = 10 g/L, Q= 500 mL/min, k3 = 10-9 m2, L= 6x10-15 m. Solid line: concentration; dotted line: velocity 65 5.15 Radially-averaged axial BSA concentrations and transmembrane velocities for an empty ECS, c30= 10 g/L, Q= 500 mL/min, k3= 10-9m2, Lp= 6x10-15m. Solid line: concentration; dotted line: velocity 65 ix 5.16 Transient radially-averaged axial protein distribution and axial velocities for a packed ECS, c30 = 10 g/L, Q = 500 mL/min, k3 = 5x10-15 m2, s = 1.0, L= 6x10'5 m. Solid line: concentration; dotted line: velocity 66 5.17 Radially-averaged BSA concentrations and transmembrane velocities for a packed ECS, c30= 1 0 g/L, Q 500 mL/min, k3= 5x10-15 m2, c = 1.0, Lp= 6x10-15m. Solid line: concentration, dotted line: velocity 66 5.18 Transient radially-averaged ECS protein concentrations for an empty ECS, including viscosity variations. Solid line: constant viscosity, dotted line: variable viscosity. c30= 20 g/L, Q = 500 mL/min, k3= 109m2, Lp= 6x10-15m 68 5.19 Transient radially-averaged ECS protein concentrations for a packed ECS, including viscosity variations. Solid line: constant viscosity, dotted line: variable viscosity. c30= 20 g/L, Q = 500 mL/min, k3 = 5x10-15 m2, E = 1.0, Lp= 6X10-15m 68 5.20 Steady-state radially-averaged BSA concentrations for c30 = 5, 15, 30, 50 and 75 g/L, Q= 500 mL/min, Lp= 6X10-15111. Solid line: packed ECS, k3 = 5x10-15m2, = 1.0; dotted line: empty ECS, k3 = 109m2 69 5.21 Time to reach steady-state as a function of initial BSA concentrations. Q — SOO mL/min, L= 6x10-15m. Solid line: packed ECS, k3 = 5x10-15 m2, E = 1.0; dotted line: empty ECS, k3 = 10-9 m2 71 5.22 Steady-state radially-averaged protein distributions at ECS conductivities of k3 = 10-9, 5x10-15, 10-16 and 10-17 m2, e = 1.0, C30 = 10 g/L, Q = 500 mL/min, L = 6x10-15 m 72 5.23 Transient radially-averaged BSA concentrations and ECS velocities for a tissue-like conductivity, c30= 10 g/L, Q= 500 mL/min, k3= 10'7m2, c = 1.0, L= 6x10-15 m. Solid line: concentration, dotted line: velocity 74 5.24 Steady-state radially-averaged protein distributions in the fluid phase at ECS porosities of e = 1.0, 0.5, 0.2 and 0.1, k3 = 10'7m2, c = 1.0, c30 = 10 g/L, Q = 500 mL/min, Lp= 6x10-15m 75 Acknowledgments I thank my advisors Jamie Piret and Bruce Bowen for the help they provided to create this thesis. I am also very grateful for the financial aid that I received from the Biotechnology group. It was a pleasure and a rewarding experience to work in the Biotechnology Lab; in particular I enjoyed the friendly atmosphere and the help received from all my co-workers. An extremely rich source of information was Anant, who also turned out to be a very good chess player. In particular I thank my "African" friend Rumina for the close friendship we developed during our graduate studies at UBC. Mitt litOen, bet. bunt 2.WciObraitg ne4cift Zott ficht Zcionerien iiiinftig ict) berMjticheit. Cured from the craving to know all, my mind Shall not henceforth be closed to pain J. W. Goethe (1825) 1 CHAPTER 1 Introduction Complex therapeutic and diagnostic proteins are the major products in mammalian cell culture biotechnology (Walker and Gingold, 1988; Atkinson and Mavituna, 1991). In conventional stirred tank bioreactors, productivity and product concentrations are low; therefore immobilized cell designs can greatly improve reactor performance (Belfort, 1989). In hollow fibre bioreactors (HFBRs), mammalian cells can grow up to 100 times higher density (>10 8 cells/mL) than in suspension type bioreactors. A typical HFBR consists of several thousand bundled ultrafiltration hollow fibres potted at each end and enclosed in a tubular shell. The cells grow around the fibres in the shell side, referred to as the extracapillary space (ECS) by analogy to the tissue-capillary system. Tharakan and Chau (1986) reported that recycle flow in the HFBR lumen creates an axial pressure drop which drives a secondary flow in the ECS referred to as the Starling flow (Figure 1.1), again by analogy to the coupled flows found in the microvascular exchange system. In the upstream half of the reactor, a fraction of the flow enters the ECS, whereas in the downstream half, it re-enters the fibre lumen (Figure 1.2). In this closed-shell, lumen recycle mode, the amount of leakage into the ECS is a function of the ECS and membrane conductivities, fibre dimensions and the pressure drop along the lumen side. Smaller molecules such as growth factors (e.g., insulin), glucose or oxygen, supplied from the lumen side, penetrate through the fibre membranes. Oxygen depletion in the ECS environment can lead to oxygen limitations and limit the viable cell distribution (Piret and Cooney, 1991). Cells are completely retained in the ECS, while high molecular proteins, introduced with the cell inoculum, are partly or completely rejected by the ultrafiltration membranes. High molecular weight products (e.g., antibodies), produced by cells and rejected by the ultrafiltration membranes are harvested from the ECS ports. The secondary 2 Figure 1.1^Pressure distribution in a closed-shell HFBR. Figure 1.2^Starling flow around a single fibre in a closed shell. 3 flow in the ECS transports these proteins towards the downstream end where they accumulate and increase the local osmotic pressure. This osmotic pressure, by reducing the total pressure difference across the membrane, can significantly affect the fluid dynamics in the ECS (Taylor et al., 1993). Transferrin mediates cellular iron-uptake and is thus an important constituent of most cell culture media. Downstream polarization of transferrin can skew cell growth towards the downstream end of the reactor which decreases overall HTER cell growth (Piret and Cooney, 1990b). On the other hand, this flow polarizes the product proteins to locally higher concentrations thus facilitating downstream processing and product purification. Recently, Taylor eta!. (1993) modelled HFBR flow dynamics and ECS protein redistribution under startup conditions where the reactor is essentially unobstructed by cells. The objective of this work was to study the transport of high molecular weight proteins in the ECS of HFBRs under packed cell conditions. To accomplish this a mathematical model was developed describing the lumen and ECS hydrodynamics in the presence of osmotically active proteins trapped in the ECS (Chapter 3). Since protein transport and hydrodynamics are coupled, numerical methods were applied. To validate the model experimentally, a packed HFBR system was designed (Chapter 4). Successful culturing of mammalian cells requires considerable technical expertise. It is also very time consuming; usually several weeks of growth are required to fill the ECS of a HFBR with cells. Therefore a packed ECS bed was simulated by filling the shell side with a gel/protein mixture (Wei and Russ, 1977). Hydrogels have a low hydraulic conductivity and are easy to prepare; upon cooling they form a homogeneous porous medium. The idea of filling the cartridge with monodisperse spherical particles (e.g., latex or glass) was rejected because such particles are expensive and it was questionable if they could be packed homogeneously between the fibres. Also, the large particle surface area could have caused significant protein adsorption. To validate the model equations, computer modelling results were compared to experimentally obtained ECS protein concentration profiles (Chapter 5). Using the model, the impact of protein loading and packed bed conductivity on the axial protein distribution and 4 hydrodynamics was investigated. By way of conclusions, the implications of the experimental and modelling results to mammalian cell culture are discussed in Chapter 6. Suggestions for further research are given in Chapter 7. 5 CHAPTER 2 Literature Review 2.1 Cell Culture In the early 1950s the first industrial use of animal cell culture was reported for the production of vaccines (e.g., polio). Since bacteria grow about 10 times faster than mammalian cells and require less expensive media, the recombinant DNA technology in the early 1970s gave hope to the idea of expressing mammalian proteins in bacteria. Unfortunately, many of these complex biological molecules require mammalian cell post-translational processing. Such protein post-translational modifications include, for example, specific modifications of amino acids or glycosylation. The development of hybridoma cell technology in 1975 by Kohler and Milstein opened the new field of monoclonal antibody production. By injecting hybridoma cells into mice, monoclonal antibodies can be harvested from the peritoneal fluid (ca. 10 mL). Due to impurities and ethical concerns, new reactor designs are desirable. Industrial production mainly uses stirred tank bioreactors, with low product concentrations and low productivities (Belfort, 1989). Additional costs arise for mammalian cell culture by the expensive media and the product recovery costs due to the presence of extraneous proteins in the solutions. Immobilized membrane reactors are attractive for culturing mammalian cells because they retain cells at low shear stress and simulate an in vivo environment. In general there are three major cell immobilization methods used, namely, attachment on microcarrier beads, entrapment in a porous matrix and physical separation by a membrane barrier (Belfort, 1989). The higher cell densities in the immobilized phase increase volumetric protein productivity and by perfusing these reactors with fresh medium the production can be maintained. After 6 startup, the cells are retained in the immobilized phase and low cell growth rates become less important. These bioreactors can be productive for weeks or even months The first use of a HFBR for mammalian cell culture was reported by Knazek et al. (1972). Cells, separated from the medium stream in the shell-side of a kidney dialysis ultrafiltration membrane module, grew to tissue-like densities and were maintained productive for nearly a month. Since then, many others have reported the use of HFBRs for culturing mammalian cells (Knazek eta!., 1974; Tharakan and Chau, 1986; Piret and Cooney, 1990a). For HiFBR the formation of necrotic regions due to nutrient limitations at high cell densities is one of the major disadvantages of cell immobilization. Also the reactor setup is more sophisticated and more training is required for its operation. Mass transfer modelling becomes an important consideration for improved reactor design, scale-up and control strategies. 2.2 Hollow Fibre Bioreactor Modelling In recent years many models have been published to describe the hydrodynamics and mass transport in HYBRs. Most of these models are based on the assumption that a multi-fibre device can be represented by a Krogh (1919) cylinder geometry where a single straight fibre, homogeneously surrounded by cells, is assumed to represent the fibre bundle. The performance of hollow fibre reactors with enzymes immobilized in the spongy section of the fibre wall was investigated by Waterland et al. (1974). A mathematical model was developed which included lumen and membrane mass transfer resistance for first-order and Michaelis-Menten kinetics. Kim and Cooney (1976) considered the same case, but reported a simpler analytical solution. Waterland et al. (1975) immobilized the enzyme I3-galactosidase in the membrane matrix and applied the previously developed model to predict the reactor performance. 7 Wei and Russ (1977) injected a gelatin/agarose gel into a HFBR shell to simulate tissue culture conductivities. A dual circuit HFBR was designed to simulate solute fluxes in capillaries. The results demonstrated that for small molecules, such as oxygen, convective fluxes can significantly contribute to their transport through tissues. Webster and Shuler (1978, 1979, 1981) considered the cases where cells or enzymes surround a single fibre and the transport occurs mainly by diffusion into the outer fibre region. Concentration gradients in the lumen, as well as the impact of the ECS hydrodynamics, were neglected. The radial diffusion equation was solved and results were presented for zero and first order kinetics in the form of effectiveness factor plots. These plots illustrated the importance of diffusive hindrances for a reacting substrate and provided a guide for the selection of HFBR dimensions. The impact of radial flow on substrate transport was first recognized by Kleinstreuer and Agarwal (1986). Cell growth as well as Michaelis-Menten kinetics of cells or enzymes in the spongy fibre matrix were considered. A model was presented that included transient lumen, membrane and ECS concentration profiles. The complete Navier-Stokes equation, governing the lumen hydrodynamics, was solved with a numerical scheme based on the SIMPLE computer code, developed by Patankar and Spalding (Patankar, 1980). However, the impact of the axial ECS velocity profiles on substrate transport was not considered. Despite the evidence of polarization from Waterland et al. (1975), who reported a 10-fold higher enzyme activity at the downstream end of a HFBR, the early models did not include convective transport effects. Tharakan and Chau (1985) measured the transmembrane pressures in a closed-shell HFBR and considered the impact of the secondary ECS flow on the cell or enzyme distributions. Models developed by Tharakan and Chau (1986) and Schonberg and Belfort (1987) only considered radial convection and neglected the influence of axial ECS transport. Libicki et al. (1988) investigated the influence of Starling flow on solute transport. They measured a 10% increase in gas transport through dense microbial pellets around a single 8 hollow fibre system enclosed in a permeable membrane. They developed a model based on the work of Apelblat et al. (1974), who modelled the fluid exchange between capillaries and tissue beds, which is similar to the flow encountered in a cell-filled ECS. An analytical solution was obtained by describing the membrane and ECS flows by Darcy's law and the lumen flow by a simplified one-dimensional version of the Navier-Stokes equation. Salmon et al. (1988) modelled the transport of reacting and non-reacting solutes for the same system investigated by Libicki eta!. (1988). Again, by incorporating the analytical velocity profiles of Apelblat et al. into the numerically solved solute transport equation, they found that, under high flow and packed cell conditions, convective transport could play a significant role. Heath et al. (1990) used nuclear magnetic resonance flow imaging to measure the ECS velocity profiles around a single hollow fibre in an empty shell. Because they were not able to place the single fibre concentrically within the cylindrical reactor shell, they numerically solved the unidirectional Navier-Stokes equation. Instead of linking the ECS and lumen flows, the maximum predicted ECS axial velocities were simply equated to the corresponding measured results. As a consequence, they obtained generally good agreement between predicted and measured velocities at other radial positions. Heath et al. (1990) also used Apelblat's solution to theoretically predict ECS velocities in both packed as well as empty multi-fibre HEBRs. However, a questionably high Darcy conductivity (10-4 m2) was arbitrarily selected for the empty ECS case. The pressure and flow distributions in a HFBR with an unobstructed ECS were modelled for different operating conditions by Bruining et al. (1989) and Kelsey et al. (1990). They assumed that radial pressure gradients in the lumen and ECS did not exist, due to the low Krogh cylinder radius to fibre length ratio. The coupled pair of second-order ordinary differential equations for the lumen and ECS pressures were solved analytically. The insulin response of a HFBR bioartifical pancreas was modelled by Pillarella and Zydney (1990) by superimposing a solute mass balance on the hydrodynamic solution of Kelsey eta!. (1990). It was shown that the radial and axial convective fluxes can dramatically 9 increase insulin response times. Their numerical modelling results agreed with transient insulin measurements. The impact of the ECS secondary flow on mammalian cell distribution and protein concentrations was investigated experimentally by Piret and Cooney (1990b). A freezesectioning technique was developed to measure the axial distributions of mammalian cells and proteins in the ECS of anisotropic membrane HFBRs following various operating protocols. They found that the downstream polarization of the high molecular weight proteins, needed for mammalian cell growth (e.g., transferrin), correlated well with the measured non-uniform axial cell growth in the HFBRs. By periodically alternating the lumen flow direction, protein polarization was reduced and a more uniform cell distribution was obtained, resulting in a 2to 3-fold increase in reactor antibody productivity. Piret and Cooney (1991) developed a model to describe oxygen gradients in HFBRs. By an order-of-magnitude analysis they showed that convective oxygen transport in the ECS can be neglected. Hence, the governing equations were simplified and an analytical solution was derived. The model included membrane and lumen resistance, as well as axial and radial oxygen gradients. The model results were presented in the form of effectiveness factors from which one can easily calculate oxygen limitations for a given reactor system. Their experimentally determined axial cell distributions agreed reasonably well with the model results. The tissue-like cell densities obtained in HFBRs can produce such high metabolic rates that gradients in nutrient concentration occur. Sardonini and DiBiasio (1992) supported their conclusions by measuring the radial mammalian cell distribution around a single hollow fibre for different lumen oxygen concentrations. None of the models described above included the impact of high molecular weight osmotically active molecules retained by the fibre membranes. Recently, Taylor et al. (1993) developed a model to describe the coupled protein transport and fluid dynamics for a Krogh cylinder geometry with an osmotically active protein solution retained in the ECS. Numerical techniques were required to solve the two-dimensional transient concentration profiles in the 10 ECS. They found that the downstream polarization of proteins can significantly change the flow distribution in HFBRs. Patkar et al. (submitted) experimentally investigated protein polarization in a HFBR and presented a simplified one-dimensional model based on the more complex two-dimensional model by Taylor et al. (1993). The transient as well as steady-state ECS protein concentrations, recovered from the frozen reactors, agreed with the model predictions. A high initial ECS protein loading reduced the non-uniform steady-state protein distributions caused by the Starling flow. Also the influence of flow cycling on the redistribution of proteins was investigated. 2.3 Protein Transport in a Gel Hydrogels are networks of organized polymers that absorb large quantities of water while remaining insoluble. These gels, which are neither solid nor liquid, have a low tendency to adsorb proteins (Gin et al., 1990). Agarose gels are made from purified linear galactosecontaining aerogel-xerogel hybrid colloids either isolated from agar or directly recovered from agar-bearing marine algae (Dean et al., 1985). The structural integrity, hardness and porosity of an agarose gel depends on the secondary structure caused by non-covalent bonds (hydrogen bonds) between various agarose chains (Dean et al.). The surface porosity of an agarose gel, observed by electron microscopy, showed that a 3% agarose gel has a very uniform internal pore structure with pore diameters ranging from 0.1 to 1.0 i_tm (Bassi et al., 1987). Therefore the agarose pore structure should offer insignificant hindrance to the transport of BSA molecules (Hjerten, 1962), which have an average diameter of about 6.3 nm (Tanford, 1961). 11 The flow field in porous media can be described by Darcy's law as long as the Reynolds numbers do not exceed a value between 1 and 10 (Scheidegger, 1960) The Reynolds number in a porous medium is defined in terms of the pore diameter of the packed bed, Re = u 'pore (2.1) where u is the superficial velocity, d is the pore diameter and u is the fluid kinematic viscosity. 2.4 Porous Bed Conductivity The presence of cells in the ECS of a HTBR can lead to significant changes in the flow behaviour and protein transport compared to an empty ECS. It is expected that, for a packed cell ECS, the protein polarization will strongly depend on the conductivity of the porous medium. Therefore, it is important to use a conductivity in the model that corresponds to the packed beds employed in the experiments. Unfortunately, conductivity values for packed beds of mammalian cells have not been reported in the literature. Libicki (1985) studied the Darcy permeability of packed spherical particles in the ECS of HFBRs and compared the results with a model by Sangani (1982). Libicki's results and the conductivities obtained from Equation 2.2 are summarized in Table 2.1. As expected, the conductivity strongly depended on the particle size and the bed porosity. Both models describe the bed conductivities reasonably, but for the highest particle radius of 87.5 Jim Sangani's modelling result agreed better with Libicki's conductivity measurement. However, it is not clear whether the conductivity of a packed bed of solid spheres is representative of the conductivity that exists for high density packed mammalian cell aggregates. Humphrey et al. (1985) reported Darcy permeabilities for bacteria in filter cakes in the range of 10"" to 12 1046m2, depending on the applied pressure difference. Due to the deformability of cells, these values are probably not valid for the low pressure gradients that exist in the ECS of I-EFBRs. Table 2.1 Conductivities of packed spheres. ^ ^ ^ Estimated Kozeny Particle Radius^Measured Porosity ^ Conductivity Conductivity Carman (Libicki) ^ (m2) (iirn) (Sangani) (m2) ^ (Eqn. 2.2) (m2) (-) 8 2.85 x1042 2.20 x 10-12 2.12 x 10-12 0.61 10 9.68 x1042 2.35 x1042 6.0 x 1043 0.44 19.8 15.9 x1042 49.0 x10-12 11.8 x 10-12 0.60 87.5 403 x10-12 447 x1042 103 x 10-12 0.52 The hydraulic resistance of packed beds of red blood cells was measured and modelled by Zydney et al. (1986). They found that these cells were compressible, therefore the packing and hence the conductivity were a function of the applied pressure gradient. A porosity E of 0.2 was measured in packed beds of cells sedimented gravitationally. The lowest porosity for uniform rhombohedral packed spheres is 26% (Bear, 1972). Based on the Kozeny-Carman equation (Kozeny, 1953; Carman, 1937), k= (2.2) 5(';(1-8))2 where S is the area of the particle surface and V is the particle volume, Zydney et al. calculated conductivities for packed, relatively uncompressed, red blood cells from 7x10-15 down to 3x10-18 m2 for the highest cell packing 13 A simple calculation of the conductivity of a mammalian cell bed based on the KozenyCarman equation was carried out here. Assuming an average mammalian cell radius of 12 1.1m, a spherical shape and packing porosities of 26, 10 and 5%, cell bed conductivities of 10-13, 4x10-15 and 4.4x10-16m2, respectively, were obtained. The flow through interstitial fibre matrices was studied by Levick (1987). He found that, even for a very highly packed cellular tissue, the flow resistance is not due to the drag effect of the bounding cell walls but from the drag created by the intracellular matrix. This intracellular matrix is created by a network of relatively coarse fixed elements, the collagen fibrils, and a finer meshwork of glycosaminoglycans. Swabb et al. (1974) investigated the hydraulic permeability for flow through subcutaneous and hepatocarcinoma rat tissue and obtained values of 6.4x10-18 and 31x10-18 m2, respectively. These conductivities are similar to the values reported by Apelblat et al. (1974) for rabbit omentum (2.6x10-18- 6.0x10-18 m2). A summary of the relevant cell bed hydraulic conductivities is given in Table 2.2. The conductivities in HTI3Rs packed with mammalian cells are expected to lie between those calculated for beds of 10 to 15 tim packed spheres and values reported for tissue regions. Table 2.2 Hydraulic bed conductivities Porous Medium^Hydraulic Conductivity^Reference (m2) Packed spheres^10-13- 4x10-16^Equation 2.2 Red blood cells^7x10-15- 3x10-18^Zydney eta!. (1986) Rat tissue^6.4x10-18, 3 1x10-18^Swabb eta!. (1974) Rabbit tissue^2.6x10-18- 6x10-18^Apelblat et al. (1974) 14 CHAPTER 3 Model Formulation Figure 3.1 shows a representative single fibre of a multi-fibre HFBR cartridge, where L is the permeable length of the fibre, Pin and P0 the pressures at the permeable fibre inlet and outlet, respectively. The fibre is potted into epoxy at both ends where no fluid exchange with the ECS occurs. The notation specifying the lumen, membrane and ECS regions as 1, 2 and 3, respectively, is used in the subsequent sections, which describe the equations governing the hydrodynamics and protein transport in all three regions. Figure 3.1 Single fibre Krogh cylinder geometry. 3.1 Krogh Cylinder Approximation To describe the flow in a reactor that consists of several thousand fibres, some major simplifications are necessary. The reactor is modelled as containing N parallel, evenly-spaced fibres that do not exchange fluid or proteins with each other. Therefore each fibre can be 15 assumed to be surrounded by an average volume of homogeneous ECS (Figure 3.1) with a hydraulic conductivity k3 and radius R3 (Figure 3.1). The radius R3, referred to as the Krogh cylinder radius, can be easily calculated from the following equation: R3 = (3.1) r- N where Rc denotes the cartridge radius. The Krogh cylinder (1919) approach was originally developed to describe the flow in analogous tissue/capillary systems, where the capillaries were assumed to be parallel and uniformly spaced. Krogh cylinder fibre lumen fibre membrane EC S/shell- side Figure 3.2 Schematic of the Krogh cylinder representation of a multi-fibre reactor. 3.2 Hydrodynamics 3.2.1 Flow in Porous Media The conductivity of a porous material defined by Darcy's law (1856) is a macroscopic property of the material. Many attempts have been made to relate the conductivity to the 16 porous structure, but no satisfactory solution has been found (Scheidegger, 1960). The Kozeny-Carman equation (Kozeny, 1953; Carman, 1937), mainly used for filtration application, models the porous medium as consisting of capillary bundles and therefore applies a classical hydrodynamic approach. The differential form of Darcy's law in vector notation can be written as k ^ V =-- •VT (3.2) where V is the superficial fluid velocity vector, t is the fluid viscosity, k is the medium conductivity and the flow potential IP is defined as follows to account for the hydrostatic pressure: qi= P+pgz,^ (3.3) where P is the pressure, g the gravitational constant and z the height with reference to a selected datum. It should be noted that, for a non-homogeneous medium, k is a second-rank symmetrical tensor. In the case of an incompressible fluid (i.e., constant p), the volume is not altered by pressure changes. Hence mass is conserved and one can write (v v), o ^ (3.4) By substitution of the isotropic form of Darcy's law into Equation 3.4, the following equation is obtained: V •(—V P), jt (3.5) Usually the isotropic conductivity k, gravity g, density p and viscosity g_t are considered to be constant. Therefore the partial differential equation (PDE), Equation 3.5, reduces to the Laplace equation, namely 17 \72 P = 0.^ (3.6) 3.2.2 Membrane Hydrodynamics The membrane is treated as a homogeneous porous cylindrical shell with a constant conductivity. Therefore the Laplace equation, which describes constant density and constant viscosity flow in any general porous medium, can be applied (Collins, 1961). Written in cylindrical coordinates and assuming axisymmetry, Equation 3.6 becomes 1 a ( a P2) a2 P2 0 + 2= r ar^ar^az (3.7) Since the aspect ratio of the membrane, (R2-R1)/L, is less than 104, pressure gradients in the axial direction are generally orders-of-magnitude smaller than in the radial direction. Thus Equation 3 7, which applies at any axial position, simplifies to 15 ( a P r^2 j= 0. r Sr ar (3.8) This ordinary differential equation is subject to the following boundary conditions: 1) P2 = (z)^(z, Ri)^at^r 2) P2 =P3(z,R2)—an3(z,R2)^at^r = R2 R, Proteins, partially rejected by the membrane, can exist in the ECS and lumen of HTBRs The boundary conditions include the osmotic pressures, III and 113, whose difference acts in the 18 opposite direction to the hydrostatic pressure driving force across the membrane. The osmotic pressure difference is multiplied by the reflection coefficient a. If a is 0, the membrane is completely permeable to the solute and solvent; if a is 1, the membrane retains the solute completely. The first boundary condition also reflects the fact that, in the simplified hydrodynamics assumed for the lumen (see Section 3.2.3), only the radially-averaged hydrostatic pressure, Pi(z), is required. By integrating this ODE twice and applying the boundary conditions, the following pressure distribution in the membrane, as a function of the local pressure difference between the lumen side and ECS, is obtained p2(r),p1(z)+ P3(Z)R2)- PI(Z)+Oirli(Z,A)-r13(2.,R2)} in 1 ^(3.9) After differentiating Equation 3.9 with respect to the radius and substituting the result into Darcy's law, the relationship for the local superficial membrane velocity becomes LR v 2(z ,r)^13,(4— P3(z , R2) +cy [113(z,R2)— II,(z,R1)] (3.10) where v represents the radial velocity and the hydraulic permeability Lp is defined (Taylor et al., 1993 ) as k2 (3 11) R, in Equation 3.10 is generally referred to as Starling's relationship, in recognition of Ernest Starling, an English physiologist who proposed in 1896 that transmembrane flow between the 19 interstitial space and capillaries can be described by a balance between hydrostatic and osmotic pressure. The resulting flow in the extracapillary space is therefore sometimes called Starling flow. For the cuprophan membranes investigated here, which have a molecular weight cut-off of about 13k Da (manufacturer's specifications), the 69k Da protein bovine serum albumin (BSA) does not leak from the ECS through the membrane into the lumen side where also no osmotic pressure exists. For these circumstances Equation 3.10 simplifies to v(z,r) =^P3(z , R2) +113(z , R2)1. 4 r (3.12) 3.2.3 Lumen Hydrodynamics The following equations of motion and continuity, written in vector notation, govern the flow in the lumen side: avi +(y•v)vi,---1 vpi+vv2v, at +Fi^ (v. v3= o^ (3.13) (3.14) where V1 represents the lumen velocity vector, I the time, p the constant density, P1 the lumen pressure, F1 the body forces and v the kinematic viscosity of the fluid. To further simplify the problem some assumptions are required. Since the fluid density is constant the body forces F1 can be neglected. Because the fluid is nearly incompressible, the laminar flow (Re < 10) in the lumen adjusts almost instantaneously to any pressure changes, and consequently, its governing equation is not time dependent. The inertial terms of the Navier Stokes equation can be 20 neglected because the radial leakage flow through the membrane is orders-of-magnitude lower than the lumen flow rate (Berman, 1953) and the flow is fully developed before it reaches the permeable region of the fibre bundle. The low fibre aspect ratio, RA, (< 10-4) allows the radial pressure gradients and the entrance effects to be neglected. In addition, the problem is axially symmetric around the fibre centre. Consequently, the radial component of Equation 3.13 simplifies to °PI — ar = (z) or while the axial component becomes OP,^1 a ( r au, (3.15) az^r ar^ar)• The last equation can be integrated twice subject to the following boundary conditions u 1^ 1= 0 at r= R1^(i.e., no-slip at the membrane wall) and a ^ Or =0 at r =0^(i.e., symmetry at the centre of the fibre) to obtain the radial velocity in the lumen as a function of the local pressure gradient, namely 2 ^u,(z,r)=—RI2 dP ( r 41=1. dz^R2) (3.16) By integrating this equation over the fibre cross section one can obtain the local radiallyaveraged axial lumen velocity as 21 1 12 rcu (zd z = R'2 d , R12 0^ 8t d z (3.17) Applying a mass balance on a differential axial control volume in the lumen and substituting Equation 3.12 into the result leads to =2 v2 (z ,R1) 2 d z^R,^R, „^, ^ PI(z)—Pjz,R2)+II3(z,R2) ). (3.18) Finally, by differentiating Equation 3.17 and substituting for the 1.h.s. of Equation 3.18, a second-order ordinary differential equation governing the lumen pressure distribution can be written, i.e., d2^ = 16LP^ ( d z2^R; (z) P3(Z,R2 )± r13(Z, R2 ) ) (3.19) This equation is subject to the following pressure boundary conditions at the fibre inlet Pin at z = 0 and fibre outlet P1= P0^at z = L One should note that the applied no-slip assumption at the membrane wall is not valid for very porous materials. In our case the membrane has a very low conductivity and the axial pressure gradients are small. Therefore the axial velocities within the membrane are very small and the axial slip at the membrane interface can be neglected (Taylor, 1971; Saffman, 1971). 22 3.2.4 Extracapillary Space Hydrodynamics For modelling purpose it is assumed, that the porous medium, such as a bed of mammalian cells, is distributed homogeneously in the ECS and therefore has a constant conductivity k3. The Reynolds numbers in a packed cell or gel ECS bed are low (Re << 1) and Darcy's law can be applied. Also, assuming axial symmetry, the problem can be described as two-dimensional and Equation 3.6 becomes I a (ra P3)-fa2 P3 = 0 r Or^ar^az2 (3.20) The fluid viscosity does not appear in this equation, because it is assumed that, over the concentration range of interest, viscosity is essentially independent of protein concentration. Equation 3.20 is subject to the following boundary conditions: 1) aP3 —0^at^z= 0, R2 r R3 2) a 133 0^at^z=L, R2 r 3) aP3 —0^at^r = R3, 0.z._. L 4) az az R3 ar aP3^R, (p, z„ „ P R2 J+II3 (z, R2 ) ) at r R2 0 or^k3 R2^1^3 The first three boundary conditions assume that the velocities at the outer boundaries of the Krogh cylinder are zero (i.e., no flow passes from one fibre space to another nor does it penetrate into the potting material at either end of the ECS). The fourth boundary condition arises from the fact that incoming and outgoing fluxes at the ECS/membrane interface must be 23 equal. After the ECS pressure distribution is obtained the axial and radial velocities, required for the protein transport equation (Section 3.3.2), can be calculated using Darcy's law. 3.3 Protein Transport 3.3.1 Protein Properties Generally the osmotic pressure for dilute low molecular weight solutions can be described by van't Hoff s (1887) law II=c RT^ (3.21) where c is the concentration of the solute, R is the gas law constant (8.314 J/mol-K), T the absolute temperature and M the molecular weight of the solute. For molecules with high molecular weights (e.g., proteins) the observed osmotic pressure does not follow this linear relationship, especially at higher concentrations. A more general form of this equation is given by a power series (Tombs and Peacocke, 1974) ^ (c+A2c2 +A3c3). ^ (3.22) The constants A2 and A3 are referred to as the virial expansion coefficients and they account for the non-ideal behaviour of the solution due to molecular interactions. The osmotic pressure relationship for BSA was extensively studied by Vilker et al. (1981a). A complex semi-empirical relationship was developed to describe the BSA osmotic pressure as a function of the protein and salt concentrations. In addition to the two virial coefficients A2 and A3, another term was included to account for the non-ideal behaviour resulting from ionic interactions, i.e., 24 RA4T c + c, A ( +,43c1+RT 2 \ Zc 1 2 N2 2 2M — 2 ins (3.23) where Z represents the BSA charge number (Z = -20.4 at pH 7.4), ms is the molar salt concentration in the solution and M is the average molecular weight of BSA (69k Da). The virial expansion coefficients are given by: A2 =-5.625x 10-4 —2.410x 10-4Z-3.664 x10-5Z2 A3= ^ 2.950 x 10' —1.051 x10-6Z— 1.762 x 10--7Z2 ^ (3.24) (3.25) Equation 3.24 and 3.25 were obtained by Vilker et. al. (1981a) by curve-fitting their measured osmotic pressure data. Vilker's et al. expressions were based on measurements at high protein concentrations (84 to 475 g/L). Since the concentration range encountered in this work was below 100 g/L, a more appropriate osmotic pressure relationship was obtained by using osmotic pressure data reported by Scatchard et al. (1946) and finding constant values for the two virial expansion coefficients by curve-fitting. In comparison with Vilker's equation, Scatchard's data show a slightly higher osmotic pressure in the 0 - 100 g/L range, a trend that was also consistent with our own measurements. In Figure 3.3, the fitted curve for Scatchard's data is shown along with Vilker's correlation. By fitting Scatchard's data to Equation 3.22, the following parameters were found which apply to a BSA concentration range of 0 - 165 g/L: A2 = 10.473 x10-3 m3/kg A3 = 17.374 x10-6 m6/kg2 25 Tombs and Peacocke (1974) also fitted experimental data to obtain a similar second-order polynomial for the BSA osmotic pressure relationship. Figure 3.3 Osmotic pressure of BSA protein solutions. Dashed line: Vilker's correlation (Equation 3.23), •: Scatchard's data, solid line: fit of Equation. 3.22 to Scatchard's data. The following simple linear relationship for the density of BSA solutions, which applies to concentrations up to 540 g/L, was found by Vilker eta!. (1981b): p = 2.54 x 10' c + 1.^ (3.26) Tanford et al. (1956) investigated the viscosity of BSA solutions for a wide pH range and expressed their results in terms of the power series = 1 + DA] c + K [pd2 c2^ (3.27a) where u is the viscosity of the BSA solution, po is the viscosity of the solvent (i.e., water at 20 °C) and [p.] is the intrinsic viscosity obtained by extrapolation of 26 11— I-to (3.27b) lim C-30 C I-t0 The constant K depends on the molecule charge, ionic strength as well as the pH of the solution. At pH 7.3, the dimensionless constant K is 2.085 and [4] = 3.65x10-3 Lig. The concentration dependent diffusivity of BSA was investigated by van den Berg and Smolders (1989). They reported that, according to other cited references, the diffusivity of BSA depends only insignificantly on the concentration of the solution. At 20 °C, the value of the diffusivity is essentially constant (D = 0.69x10-1° m2/s) over a wide concentration range (0 to 300 g/L). Anderson et al. (1978) investigated the diffusivity at pH 6.5 and obtained by extrapolation the following linear relationship (3.28) D = D0(1+ 6x10 -4 c)^ where the Do is the diffusivity at infinite dilution (0.59x10-1° m2/s at pH 6.5) and c has units of g. In the range of concentrations of interest in this study (up to 100 g/L), neither density nor diffusivity are strongly concentration dependent (Table 3.1) However, the viscosity does show a significant concentration dependence and the impact of this will be discussed further in a later section. Table 3.1 Influence of BSA concentration on viscosity, density and diffusivity; percent changes in value from an infinitely diluted solution CBSA^ 50 WI' CBsA = 100 g/L density^p^(Eqn. 3.26) 1.3% 2.6% diffusivity^D (Eqn. 3.28) 3% 6% viscosity^IA^(Eqn. 3.27a) 25% 64% 27 3.3.2 Governing Equations For the axially symmetric ECS region, the equation of protein continuity for a constant density p in the absence of any sinks or sources can be written as follows: 1 5c 3 = -a--(K Drac3)+ a (K Dac3) K 113 ac3 -K ^ az^aZ^c^OZ^c^ar at r^d^d (3.29) where c represents the protein concentration in the fluid phase, c the packed bed porosity, u the axial and v the radial fluid velocity, D the protein diffusivity and t the time. For completeness, the two protein hindrance factors Kd and IC for diffusive and convective transport respectively, are included. Since the BSA molecule is about two orders-ofmagnitude smaller than the average size of the gel pores, no hindrance is expected and a value of 1 for Kd and IC was assumed. The partial differential equation is subject to the following boundary conditions: 1) c3-c30^at^0_<_z.L,^1=0 2) ac 3^ =0 az at^z =0,^R2 .t- R3, 1 0 3) ac3 0 az at^z= L,^R2 _r_ R3, 1 0 4) 3 =O ar at^0< z < L,^r R3,^> v c3-Kd^at^0< z < L,^r =R2,^t> O . 5)^Kc-l 6^ar The first boundary condition assumes that initially the proteins are distributed homogeneously throughout the entire ECS. The subsequent conditions state that there is no protein leakage flux through any of the ECS boundaries. Boundary condition 5 implies that, at the membrane 28 surface, since there is no protein leakage, the convective and the diffusive protein fluxes must be in opposite directions and equal in magnitude. 3.3.3 Simplified One-Dimensional Model In this section a simplified one-dimensional (i.e., radially-averaged) version of the model is developed to describe the hydrodynamics and protein transport in the packed ECS surrounding a single fibre. It is expected that, because of the typically low aspect ratios ((R3-R2)/L < 10-3-10-4) of the ECS, such a one-dimensional model should reasonably approximate the fully two-dimensional model developed above. Applying a mass balance to a differential axial control volume in the ECS results in dri3^2R3 ^ vAz R,,) d z (R; — R22 ) (3.30) where the local axial velocity (ii3) is averaged over the radial cross-section of the ECS By differentiation of Darcy's law, a second relationship for the velocity gradient results, namely d113^k3 d2 P3 (3.31) d z^j.idz2 The local radial velocity at the junction of the ECS and the fibre wall can be obtained from Equation 3.12. Substituting Equations 3.12 (with r = R2) and 3.31 into Equation 3.30 yields the following governing equation for the pressure distribution in the ECS: d2 p3^2R1^L dz2^(R.:—R)k3 (pi_ p3+113) (3.32) 29 Because of the closed shell operating condition, no fluid can leave the boundaries at z = 0 and L, and one can write: d1 ,00 dz d P3 0 dz at^z = 0 at^z = L . In the one-dimensional model for an empty ECS, the conductivity k3 is given explicitly by the following expression (Taylor eta!., 1993): k3 = 1( 4R 34^R log --J-)-3R32 +R.; 8 1Z —^\R2 (3.33) In the absence of osmotically active molecules, the coupled second-order differential equations (Equation 3.19 and 3.32), governing the lumen and ECS pressures, can be solved analytically. The solution is presented in Appendix 2. The one-dimensional protein transport in the ECS can be described by the simplified diffusion-convection equation a c3 a (ICD at^OZ^ where U3 - a (K azal^az c 6 3 (3.34) represents the local radially-averaged ECS protein concentration. The radially- averaged ECS axial velocity 173 is not constant and therefore must remain inside the differential. Equation 3.34 is subject to the following boundary conditions: 1)^Z73 = C30^ at^0 < z < L,^/=0 30 2) ac3 =0 at^z = 0,^t 0 3) ac3 - 0 at^z = L,^t az az The first boundary condition reflects the initial homogeneous protein distribution in the ECS The other two boundary conditions state that there is no protein transport through either end of the ECS. 3.4 Numerical Methods The coupled hydrodynamic and protein transport equations were solved numerically by finite difference methods. The spatial derivatives in all equations were discretized using the control volume method recommended by Patankar (1980). A listing of the FORTRAN source-code for the one-dimensional and two-dimensional computer models is included in Appendix 4. The system parameters were read from an external input file, all variables were specified and the program was started with a uniform initial protein and osmotic pressure distribution. The partial differential equation that governs the pressure distribution in the ECS was solved by a line-by-line scheme with over-relaxation (Patankar, 1980) to ensure convergence. After each "sweep" of the ECS pressures in the axial and radial directions, the lumen pressures were updated until the maximum pressure change between two complete "sweeps" was less that 10-12 Pa. Due to the extremely low membrane leakage flow, the lumen pressure did not change significantly. It was subsequently found, therefore, that the lumen pressures only needed to be updated after the ECS pressures had converged. After calculating the ECS pressures, the radial and axial ECS velocities were calculated halfway between the pressure nodes using a central difference approximation of Darcy's law. 31 The alternating direction implicit (ADI) scheme developed by Peaceman and Rachford (1955) was used to solve the discretized transient form of the two-dimensional protein diffusion-convection equation. This scheme is the recommended one, if applicable, because it requires minimal computational effort. Locally high Peclet numbers can create tridiagonal matrices which are no longer diagonally dominant. Therefore, the first-order upwind corrected "power law" scheme proposed by Patankar (1980) was included. In the simplified one-dimensional version of the program, the protein diffusion-convection equation was solved by a Crank-Nicolson method (Patankar, 1980) with a similar upwinding scheme. Since the osmotic pressure changes associated with the new protein distribution in the ECS were relatively small, the velocities were at a quasi-steady-state and, consequently, the velocity field was lagged one time increment behind the concentration field. These assumptions hold only as long as the flow field changes over each new time step are relatively insignificant. The osmotic pressures, associated with the local protein concentrations at the membrane/ECS boundary, were updated and the new hydrodynamic conditions calculated. At pre-defined time steps all information was written to separate files until the maximum time was reached or the slope of the maximum concentration difference between two successive time steps was less than 10-8g/L-s. To check for protein losses, associated with numerical instabilities, a mass balance was occasionally calculated by integrating over the permeable length of the fibre a set of third-order splines fitted to the radially-averaged ECS concentrations. The Dirichlet boundary conditions associated with the lumen pressures did not require any special treatment to close the tridiagonal matrix. To handle the Neumann boundary conditions for the ECS pressure as well as the diffusion-convection equation, integration over half control volumes was necessary (Patankar, 1980). The equations were discretized to allow for a variable grid spacing. In Figure 3.4 a control volume is sketched showing the notation used. The control volume faces were 32 located halfway between neighbouring grid points. The indices i and j define the axial and radial grid locations and the associated grid spacing. The velocities were calculated at the control volume boundaries (dashed line); whereas the pressures and the concentrations were calculated at the nodal positions. Figure 3.4^Computational control volume. 33 CHAPTER 4 Materials and Methods 4.1 Hollow Fibre Bioreactor Specifications A gambro® GFE 15 (Hechingen, Germany) cartridge was split in half along its axis and geometric measurements were taken with a dial caliper (Mitutoyo, Japan). The empty cartridge shell volume was calculated based on these measurements. A method described by Klein et al. (1977) was used to measure the fibre volume increase during wetting. A light pharmaceutical mineral oil (Stanley Pharmaceutical, Vancouver BC) was injected into the cartridge shell. Since the hydrophilic fibre material does not extensively adsorb the hydrophobic oil, the unwetted cartridge ECS volume can be calculated based on the measured oil volume (Klein et al., 1977). A water flow was then established in the fibre lumina and oil ejected from the ECS ports, due to fibre swelling, was collected and weighed. Based on the known density of oil (p = 0.852 kg/L), the ECS volume in the presence of wet fibres was then calculated. Since the ejection of oil from the ECS port occurred almost instantaneously and also no further oil flow was noticed up to about 5 min afterwards, it was assumed that water leakage into the ECS was minimal during this time. 4.1.1 Fibre Dimensions The hollow fibre ultrafiltration membranes in the gambro® GFE 15 reactor consist of cuprophan, a cellulose material made from cotton linters manufactured via a cuprammonium intermediate (Klein et al., 1977) Unlike hydrophobic polyamide or polysulphone membranes this hydrophilic material adsorbs little protein. Cellulose based membrane materials change dimensions when wetted with water, which is referred to as swelling (Klein et al., 1977; 34 Broek et al., 1992,). Therefore dry and wet fibre dimensions can differ significantly. By placing dry and wet fibres under a light microscope (Leitz; Wetzlar, Germany), the outer fibre diameter was measured in each case. It was not practical to cut a single wet fibre to determine the inner radius. Since the fibre membrane permeability was extremely low, the lumen inner diameter was therefore estimated based on pressure drop measurements in the HFBR cartridge using the Hagen-Poiseuille relationship for laminar flow in tubes —P nR4 Pin^out ^ 8p.^L (4.1) where Q1 is the volumetric flow rate through a single fibre. The pressures at the upstream and downstream lumen manifolds were measured by inserting small syringe needles through the cartridge shell and using water manometers. Needles were positioned at several locations near the potting surface of the manifold to measure the radial manifold pressure gradients. The pressure drop in the fibre potting region was measured using a cartridge end sawed at the fibre bundle side of the potting surface. The pressure at the manifold side of the potting was determined for several volumetric flow rates. For the latter measurements, the pressure at the potting outlet was assumed to be atmospheric. 4.1.2 Membrane Permeability Measurements The hydraulic permeability of the GFE 15 fibre membranes was determined using a cartridge which was wetted over night. One lumen port was sealed and a constant water pressure was applied to the other port. The volumetric flow rate was measured through the ECS port closest to the sealed lumen port. It was experimentally validated that the volumetric flow rate does not depend on the location of the exit port by also measuring the flow with both ECS 35 ports open. The hydraulic permeability Lp was then calculated from the following relationship (Mulder, 1991): J= Q L = AP Amem^p, (4.2) where J denotes the membrane flux, Q3 the volumetric flow rate, Amen, the inner membrane surface area, p. the fluid viscosity and AP the applied pressure difference between the lumen inlet and ECS outlet. The membrane surface area was calculated from the measured fibre dimensions and was in close agreement with the value supplied by the manufacturer. When using Equation 4.2, it was assumed that all of the measured pressure drop occurred across the membrane. This assumption was reasonable since the measured volumetric ECS flow rates were extremely low (maximum 6 mL/min) and hence both the lumen and ECS axial pressure losses were negligible. 4.2 Gel Preparation To avoid protein heat denaturation, a low melting temperature (ca. 30 °C) agarose gel type VII (Sigma, St. Louis MO) was chosen to fill the ECS of the reactor. The gel showed sufficient strength at 2 - 3% (w/w) to retain its integrity. To prepare the gel, PBS buffer (pH 7.4) was heated to 60 °C before a weighed amount of agarose powder was added. A clear viscous liquid was obtained after 1 - 2 h of stirring. A solution of bovine serum albumin (BSA) (Sigma) was prepared in PBS with the addition of 0.2% sodium azide to prevent bacterial growth. Trace amounts of azoalbumin (Sigma) were added to the BSA solution in some experiments. Azoalbumin, a BSA derivative with azo-groups, gives the protein solution a light reddish colour at neutral pH. Equal amounts of the protein and agarose solutions were 36 mixed together after the gel had cooled to about 40 °C. After preparation the final gel/protein mixture had half of the initial BSA and agarose concentrations. 4.3 Reactor Preparation and Flow Configuration The lumen, as well as the ECS ports, of a dry HFBR cartridge were closed with rubber tubing before the reactor was placed into a 37 °C warm room. After about 2 h the cartridge was at equilibrium with the warm room temperature and ready for filling. The previously prepared protein/gel mixture was then added with a peristaltic pump (Amicon, Danvers MA) from the bottom ECS port into the HFBR which was tilted 45° from vertical. A very slow volumetric flow rate (10 mL/min) was used to avoid excessive formation of air bubbles in the gel. After the ECS was completely filled with gel, the ECS flow was recycled at 50 mL/min until most of the remaining small air bubbles were removed and a homogeneous gel obtained. The reactor was then taken out of the warm room for gelling at room temperature (20 °C). The following day the ECS tubing ports were carefully cut off with a hack-saw and afterwards the two openings were sealed with a thin layer of water resistant epoxy glue to prevent leakage. Figure 4.1 shows the HFBR flow configuration for the protein polarization experiments. A constant head of PBS buffer (with 0.2% of sodium azide) was established. After the liquid passed through the lumen side of the horizontal reactor, it was returned to the head vessel with a peristaltic pump (Amicon, Danvers MA) from a vertical tube in the recycle vessel. This configuration had the advantage of creating a constant non-pulsating lumen flow. Parafilm covers (American National Can, Greenwitch CT) on the head and recycle vessel prevented water evaporation during the long runs. The total lumen volume in each experiment was about 1 L of buffer solution. Liquid volumes, collected at the reactor outlet, were measured to determine lumen flow rates. Only short collection times (maximum 20 s) were used to 37 avoid substantial changes in lumen pressures. Afterwards the collected liquid was poured back into the lumen recycle loop Each measurement was repeated at least twice. After a defined polarization time the flow was stopped and the lumen fluid replaced with air before the whole reactor was frozen in liquid nitrogen. The advantages of freezing were that the fragile gel in the reactor solidified and protein redistribution during cutting was prevented. Figure 4.1 Schematic diagram of the reactor flow configuration. 4.3.1 Preparation and Analysis of Samples After the reactor was cut with a hack-saw into several (5 - 11) sections, the frozen gel was removed from the shell, weighed and diluted with a defined amount of water (about 10-fold). The thawed gel and fibres were broken to decrease the time required to obtain a uniform protein concentration in the gel/water mixture, since the protein transport from the gel occurred mainly by diffusion. The gel/water suspensions were occasionally stirred and 38 samples were taken. After about 2 days, the protein concentrations no longer varied, indicating that an equilibrium had been reached between the protein concentrations in the suspended gel particles and the water. Also, initial protein concentrations, based on weight measurements, agreed with these measurements. To measure the protein concentrations in the aqueous solutions, a small volume of liquid (1 mL) was taken with a pipette such that no fibres or gel particles were sampled. The sample was spun for 15 min in an Eppendorf centrifuge (Eppendorf 5415; Hamburg, Germany) at 14,000 rpm (ca. 800 g) to separate small agarose or fibre particles that could influence the protein measurements. A Bio-Rad protein assay (Bio-Rad Laboratories, Richmond ON) was used to measure the protein concentration of the samples. The Bio-Rad assay dye binds to the BSA molecules and creates a blue solution. The Bio-Rad supplied BSA protein standard was used to prepare 6 different solutions of known concentration (50 to 260 mg/L). The supernatants of the spun samples were diluted with water such that the protein concentrations always remained in the range of the prepared BSA standards. After pipetting 20 IAL of the protein sample into a multi-well plate, 180 tL of five times diluted Bio-Rad dye was added. Each sample and protein standard was prepared at least in triplicate. After about 20 min the absorbances were analyzed with a Molecular Devices Vnia, microplate reader (Menlo Park, CA) at a wavelength of 595 nm. The integrated software performed a regression of the BSA standard and calculated the sample protein concentrations as well as the standard deviations. To ensure complete protein binding with the dye, the plates were analyzed again after 10 min; no further absorbance changes were measured. The reactor ECS protein concentration in the cut section was calculated based on the gel weight, assuming that the gel and the protein solution had the same density as water. Since the membrane excluded BSA, the water in the membrane as well as the dry fibre weight was subtracted (1.47 g/cm-cut) from the weight of the cut section. A sample protein concentration calculation is included in Appendix 3. 39 4.4 Gel Conductivity Measurements An agarose/water gel was prepared according to the method in Section 4.2. After the solution cooled down to about 33 °C the liquid gel was poured into a 60 mL Buchner funnel (Pyrex #36060C, Corning NY) with a coarse glass fit and an inner diameter of 4 cm. This procedure was performed carefully in order to prevent the creation of entrapped air bubbles in the gel. After the gel had solidified (3 - 4 h), the gel conductivity was measured using the falling head technique (Bear, 1988). A schematic of the this technique is shown in Figure 4.2. Figure 4.2 Schematic diagram of the falling head method for gel conductivity measurements. A constant water head was applied at the gel side while on the other side, the water level in the capillary was measured as a function of time. The gel conductivity k was then calculated based on the following equation (Bear, 1988): t /1„,^( k ^ In ^ k Aga p g (4.1) 40 where g is the gravitational constant, p is the density of water, 5 is the gel thickness, vi is the fluid viscosity, Age, and A cap are the cross-sectional areas of the gel and the capillary, respectively, h0 is the height of the constant water head, and h is the water height in the capillary at time 1. By plotting the modified water height ( ln[h .1(h 0-h)] ) versus time the gel conductivity can be calculated from the slope of the regressed straight line. To ensure a constant h., the water reservoir surface area was orders-of-magnitude greater than the capillary surface area. Since the frit permeability was also orders-of-magnitude (>106) higher than the values measured for the gels, the pressure drop in the frit did not influence the results and was therefore neglected. A volumetric pipette (0.1 mL) was used as the capillary. A ruler was fixed next to the capillary and the water heights h was measured as a function of time. To avoid gel compression due to the applied pressure, water heads of less than 20 cm were used; in this range a constant value of k was obtained independent of he at each experiment. To ensure that the water did not leak around the gel, a blue ink (Pelikan 4001, Germany) was injected into the water near the gel surface to allow visual observation of any bypassing. 41 CHAPTER 5 Results and Discussion 5.1 Reactor Specifications The geometric specifications of the reactor, the dimensions of the fibres, the number of fibres in a cartridge and the hydraulic permeability of the fibre membranes are amongst the reactor characteristics required for modelling. A total number of 8128 hollow fibres were counted manually in a gambro® GFE cartridge sacrificed for this purpose. A uniform fibre length of about 10 cm was cut from the fibre bundle of another cartridge and weighed. Fifty fibres were manually separated from this bundle and also weighed. From these two weight measurements, the total number of fibres in this bundle was estimated to be 7950. This result differs only slightly (2.2%) from, and hence confirms, the number of fibres obtained by counting. An average membrane hydraulic permeability of 6x10-15 m (± 10%) was determined from measurements made on several cartridges. The variations were considered to be insignificant and it was assumed that they were mainly due to changes in the membrane material rather than to varying numbers of fibres in cartridges. The manufacturer reported a slightly higher membrane hydraulic permeability of 8x10-15 m. To investigate the impact of membrane protein adsorption on the membrane hydraulic permeability, the ECS of a HFBR was loaded with BSA (10 g/L) overnight. After flushing the ECS with water on the following day to remove the BSA, no changes in permeability were measured. Solutions of BSA (1.1 and 2.2 g/L), containing different amounts of cut fibres, were prepared and the protein concentrations were measured over several days. Since no concentration changes were found, it was verified that BSA did not adsorb to the hydrophilic membrane structure. These results are very different from those found by Patkar et al. 42 (submitted) for hydrophobic polysulfone fibres, where significant BSA adsorption occurred in a few hours. The pressure drop in the fibre potting was determined separately using a sawed HFBR cartridge consisting only of the inlet manifold and the upstream potting section of the fibre bundle. From pressure and volumetric flow rate measurements in this specially prepared cartridge, an average inner fibre radius of 88 1.11T1 was calculated by using the Hagen-Poiseuille equation. An average inner fibre diameter of 91 ±1 pm in the potting material was measured under the light microscope by Patkar et al. (submitted). The small difference between these two measurements may have been due to slight swelling of these outwardly constrained fibres upon wetting or the fact that the fibre cross-section in the potting region was not perfectly circular (Patkar et al., submitted). The outer radii of dry hollow fibres were determined under a light microscope. The visually distinguishable inner membrane boundaries also allowed measurements of inner fibre radii. After wetting with water, the fibres became transparent and the inner fibre radii were no longer visible under the light microscope. An average inner fibre radius of 115 pm was calculated from the lumen pressure drop and flow rate measurements using the HagenPoiseuille relationship. Since the fibres in the inlet and outlet potting regions had a smaller average inner fibre radius, the measured pressure losses in these two regions were subtracted from the overall pressure drop. The calculated average inner lumen radius agreed well with the assumption that, during wetting, the inner and outer membrane radius increased about 15 whereas the thickness of the membrane hardly changed (ca. 1 pm). When single fibres were wetted, an average 10% increase in axial length was measured. A cartridge volume of 191 mL (without fibres) was calculated from the measurements shown in Figure 5.1. Dry and wet ECS volumes of 127 mL and 101 mL, respectively, were calculated by subtracting the total fibre volume from the HFBR shell volume. This results were in excellent agreement with the calculations made from the oil measurements for the dry *II cr.z1. C -I el) :J1 )—, 4 CD lw v) Z .-1 co ca.. a. §. ct) cn . 5 r, o CD CM Po 0 0-i 0 e C) ,-1-1 tal — Lt, CZ 4 c) A) -1. ca. Cm? co 4=, L..) 44 and wet reactor (128 and 102 mL). A summary of the measured and manufacturer supplied values is given in Table 5.1. Table 5.1 Summary of gambro® GFE 15 reactor specifications. Parameter Dry Wet Inner fibre radius 100" plin 110*, 115 Outer fibre radius 108*/** ptm 126 *, 124** pm Length of fibre 21.5** CIT1 23.6** C111 Length of potting material 1.2** CII1 - Number of fibres 8128** Inner membrane area 1.4*,^1.5** m2 Membrane permeability 8^6** x10-15 m 102** mL Total cartridge volume 191** mL ECS volume 128** mL Membrane volume 10.0*/** nth 11.8*,^13.0** nth Lumen volume 76.4*/** mL 95.7*, 93.7** mL ** manufacturer's specification measured Park and Chang (1986) investigated the radial lumen flow distribution due to radial pressure losses in the inlet and outlet manifold of multi-fibre cartridges. They reported that, under conditions where manifold pressure losses are significant, the highest lumen velocity occurred at the centre of the fibre bundle and decreased towards the periphery. A relatively 45 uniform lumen flow distribution was obtained for large axial pressure drop parameters (K), small Reynolds numbers (Re) and large axial manifold thicknesses, where K and Re are special parameters defined by Park and Chang. For a gambro® reactor, having a typical lumen recycle flow rate of 500 mL/min and an axial manifold thickness of 1.2 cm (Figure 5.1), an axial pressure drop parameter (K) of 2x105 and a Reynolds (Re) of 170 was calculated. For these conditions, a radially uniform fibre velocity distribution was expected (Park and Chang). This was verified by measuring the radial pressure distribution in the entrance manifold of a HFBR cartridge. For inlet pressures of 4000 Pa, pressure variations of less than 1% were found between the reactor symmetry axis and the outer manifold radius. 5.2 Gel Conductivity Measurements The falling head method (Bear, 1972) was used to determine conductivities for agarose gel concentrations of 1, 2 and 3%. All experiments were performed at room temperature and low pressure gradients (< 2000 Pa) to avoid gel compression. In Table 5.2 a summary of the calculated agarose gel conductivities, based on the measurements, is given. A correlation coefficient r> 0.99 for the linear 1n[(h-170)1h0] versus time plots was obtained for all conductivity calculations. The low Reynolds number flow in the porous media was described by Darcy's law. The porous medium Reynolds number (Re << 1), calculated from Equation 2.1 using the maximum superficial velocity encountered, was orders-of-magnitude below the upper limit reported for the Darcy flow regime (Scheidegger, 1960; Bear, 1972). The two conductivity measurements for a 2% gel were not identical. This might be due to a slightly altered gel structure because of higher initial pressures applied in the first 2% gel experiments (Table 5.2), or to minor leakage around the edges of the gel plug. The conductivity of the 1% gel was about 10-fold greater than that of the 2% gels, due to the lower agarose concentration. But the increase in conductivity was not as substantial when the 46 gel concentration was reduced from 3% to 2% A 2% agarose gel with a 2% BSA addition was prepared to investigate the effect of the presence of protein on the conductivity of the gel used in the gel/protein-filled HFBR experiments. It was noticed that the BSA containing gel mixture had a softer consistency than a pure agarose gel. The presence of BSA increased the conductivity of the gel, presumably because this contaminant adversely affected the gel structure. Table 5.2 Experimentally determined agarose gel conductivities. Agarose Range of Gel Gel Concentration Applied Pressure Thickness % P (Pa) 5 (cm) 3% 1577 - 1322 1.95 0.9x10-16 2% 1555 - 180 2.15 2.7x10-16 2% 930 - 852 2.05 1.5x10-16 1% 1822 - 627 2.50 20x10-16 2% + 2% BSA 940 - 700 2.30 7.5x10-16 (w/w) Calculated Gel Conductivity k (m2) The measured results did not agree with the 2% agar-agar gel conductivity of 4.4x1043 m2 reported by Pallmann and Deuel (1945). They measured the conductivity by adding the gel solution to a tube with a frit, allowing it to solidify and collecting the fluid, transmitted through the gel under an applied pressure (1 m water), for several days. Obviously, their reported conductivity values were not corrected for viscosity, or the wrong units were presented. If these values were multiplied by the viscosity of water at 20 °C, the 47 conductivities would agree with the ones found here for 2% agarose gels. Wei and Russ (1977) measured a conductivity of 2x10-'7 m2 for gelatin/agarose (3.5%/1.25%) gels. They determined the conductivity of the gel by measuring the movement of the percolating fluid in a capillary using a light microscope. This lower conductivity, given the higher overall gel concentration, is in reasonable agreement with the results measured in this study. 5.3 Model Parameters The Krogh cylinder radius R3 was calculated assuming a homogenous fibre distribution in the cartridge with a constant inner shell diameter of 3.15 cm (Figure 5.1). Since the cylindrical cartridge shell diameter increased radially to 4.12 cm over a length of 2.1 cm at both ends, the outer part of the fibre bundle was exposed to an additional ECS volume of 11.9 mL or 9% of the total ECS volume at each end. The effect of these additional spaces on ECS hydrodynamics and protein redistribution cannot be properly accounted for by the Krogh cylinder assumption and therefore it was neglected. It also was assumed that infiltration of agarose into the ultrafiltration membrane did not alter the hydraulic permeability of the fibre membranes. This is a reasonable assumption since the agarose gel conductivity was still three orders-of-magnitude greater than the conductivity of the cuprophan membrane material (5x10-19 m2). The low Reynolds number flow in the porous media was described by Darcy's law. The highest axial ECS velocity encountered during computer simulation, along with an assumed agarose gel pore diameter of 1.0 pun, yielded a maximum Reynolds numbers of 10-5, which was still orders-of-magnitude lower than the upper limit (Re < 10) reported for the Darcy flow regime (Scheidegger, 1960; Bear, 1972). The stretching of the fibres in the cartridge shell, due to 10% axial expansion during wetting, was neglected. Only a small fraction of the lumen flow passed into the ECS (<<1%) and thus a nearly linear axial pressure distribution existed in the fibre lumen. Since the 48 volumetric lumen flow rate was measured, the pressure boundary condition at the permeable fibre inlet was calculated by using the Hagen-Poiseuille relationship for laminar flow in a straight tube. The pressure at the permeable fibre outlet was assumed to be atmospheric. The fibre length was taken to be the axial distance between the two inner surfaces of the inlet and outlet potting sections (21.5 cm). Pressure losses in the impermeable fibre pottings were not considered, since these did not contribute to any ECS leakage flow. The reported lumen flow rate is, by definition, the total flow into a multi-fibre HFBR. The volumetric flow rate for a single fibre Q1 was obtained by dividing the total flow through the HFBR Q, uniformly amongst all 8128 fibres in the cartridge. An ECS porosity of 100 % (6 = 1) was used for nearly all computer modelling. Only at the last Figure (5.24) the influence of the decreased ECS porosities on the protein transport was investigated for a very low ECS conductivity. All computer simulation results presented were obtained using 201 uniformly spaced grid points in the axial direction and, for the twodimensional simulations, 11 uniformly spaced radial grid points at each axial node. In the onedimensional empty ECS model, extremely steep concentration gradients were encountered as steady-state was approached. However, 201 axial nodes were still sufficient to handle these steep gradients. Computer simulations using 801 axial nodes increased the computational time but did not alter the results. Since extreme concentration gradients were not found in the packed ECS simulations, 201 axial grid points were more than sufficient. The same steadystate solution was obtained for a wide range of time-increments; 360 - 18000 s for an ECS conductivity of 5x1045 m2 and 10 - 500 sec for an empty ECS. However, for the transient model simulations, sufficiently small time-steps were always used to ensure accurate results, 5. 1000 s for a packed and s for an empty ECS. The ECS concentration profiles were considered to be at steady-state when the slope of the maximum concentration change between two successive time steps was less than 10-8g/L-s. A summary of the computer modelling parameters is given in Table 5.3. 49 Table 5.3 Modelling parameters for HTER protein polarization program. Parameter Values Units Symbol Temperature 293 K T Permeability of membrane 6x10-15 m LP Viscosity of water 1.002x10-3 Pa-s i_t Fibre length 0.215 m L Inner fibre radius 115 iim R1 Outer fibre radius 124 1-tal R2 Krogh radius 175 jAM R3 ECS conductivity 10-9- 5x10-17 m2 k3 ECS porosity 1.0 Molecular weight of BSA 69,000 DA M Initial concentrations of BSA 5 - 75 gli- c 30 2nd virial coefficient of BSA 10.473 x 10-3 L/g A2 3rd virial coefficient of BSA 17.374 x 10-6 L2/g2 A3 Diffusivity of BSA 0.69x10-1° m2is D HTER flow rate 447 - 544 mL/min 0 Number of fibres in HFBR 8128 - N 6 50 5.4 Modelling of Hollow Fibre Bioreactor Hydrodynamics Apelblat et al. (1974) modelled the fluid dynamics in a tissue-capillary system using the Krogh cylinder approach. Their system corresponds to the packed ECS fliFBR system in the absence of osmotically-active proteins. Apelblat et al. solved the Laplace equations governing the pressure in the ECS, lumen and membrane by power series expansions. This publication, as well as others which cited and attempted to correct it (Heath et al., 1990; Libicki, 1986), contained equations with typographical errors. A corrected solution of the Apelblat hydrodynamic model is included as Appendix 1. Figure 5.2 Axial ECS pressures at R2 calculated using Apelblat's analytical solution (-) and the 2-D (.) and 1-D (0) computer models. k3= 5x10-15 M2, c= 1.0, API = 5000 Pa, L = 6x10-15 M. A program was written to calculate the Apelblat et al. solution, in terms of Bessel series, for the axial pressure distribution along the outer radius of the fibre membrane (R2) for an ECS conductivity of 5x10-15 m2. These results as well as the ones obtained from the two- 51 dimensional and one-dimensional computer models are plotted in Figure 5.2. However, Figure 5.2 reveals that pressure variations between all three model predictions were negligible at all axial positions. Thus, the one-dimensional model described the HiFBR hydrodynamics as well as the two-dimensional model and, therefore, was used to generate the results shown in Figures 5.3 to 5.5. Axial Distance (cm) Figure 5.3 Axial ECS pressure distribution for different conductivities k3, 0 = 500 mL/min, e = 1.0, Lp= 6X10-15 M. In Figure 5.3 axial lumen and ECS pressure distributions for different ECS conductivities are shown. Since the leakage flows from the lumen into the ECS were less than 0.01% of the lumen inlet flows, even in the empty ECS case, a nearly constant lumen pressure drop was encountered in all computer simulations. In the first half of the fibre, the ECS pressures exceeded the lumen pressures, resulting in a positive leakage flow into the ECS (Figure 5.4). Fluid re-entry into the fibre lumen results in a symmetric leakage flow downstream. 52 Figure 5.4 Transmembrane and average ECS velocities for different conductivities k3. Q = 500 mL/min, c = 1.0, Lp= 6x10'5 m. Solid line: transmembrane velocity, dotted line: radially-averaged axial ECS velocity. A mirror image of the pressures exists at half the fibre length where all pressures intersect. At this point no transmembrane flux between lumen and ECS occurs and the axial velocities reach their maxima (Figure 5.4). The transmembrane velocities are governed by the transmembrane pressure differences which are influenced by the ECS conductivity, membrane permeability and the lumen pressure drop. For an ECS conductivity of 10-9 m2, corresponding to an empty ECS situation, the ECS pressure distribution was nearly constant. A reduction of the ECS conductivity increased the flow resistance there and hence, the ECS pressures approached the lumen values (Figure 5.3). For the lowest plotted ECS conductivity of 5x10-15 m2, transmembrane pressure differences existed only over very narrow regions near the fibre inlet and outlet. This resulted in a nearly uniform axial ECS velocity profile (Figure 5.4). The average axial velocities for an ECS conductivity of 10-17 m2 were not plotted in Figure 5.4, since these velocities were extremely low (<0.5 tim/h). A further conductivity 53 decrease resulted in such low transmembrane pressure differences that, round-off errors, even with double precision computing, created numerical problems. Axial Distance (cm) Figure 5.5 Axial ECS pressure distribution for different membrane hydraulic permeabilities. Q = 500 mL/min, k3 = 5x10-15 m2, E = 1.0. The lumen and ECS pressure distributions for an ECS conductivity of 5x10-15 m2 are plotted as a function of membrane permeability in Figure 5.5. Decreasing membrane permeabilities resulted in increased transmembrane pressures. The flow resistance in the membrane increased and hence the resistance in the ECS became less important. At a membrane hydraulic permeability of 6x10-15 m, the main flow resistance was found in the ECS. Decreased membrane permeabilities resulted in more uniform ECS pressure profiles, like those found in the empty ECS case (Figure 5.3). An ECS pressure profile that followed the lumen pressure distribution was found for large membrane hydraulic permeability values (6x10 13 m). Here, the flow resistance in the membrane became insignificant and compared - to that in the ECS. 54 The current model neglects axial pressure gradients in the membrane region which could influence the results particularly at high membrane permeabilities or low ECS conductivities when the membrane resistance is relatively small. But, since the membrane volume is more than an order-of-magnitude less than the ECS volume, the influence of the axial membrane leakage flow on the Starling flow and protein transport in the ECS was considered to be minor. 5.5 Protein Transport in a Packed Hollow Fibre Bioreactor ECS convective flow transports initially uniformly distributed proteins towards the downstream end of HYBRs. The proteins accumulate until the osmotic pressures, associated with the local protein concentrations, match the lumen pressures. A steady-state situation is achieved after all proteins are swept downstream and distributed such that the effective (i.e., P3 - n3) ECS pressures match the lumen pressures. In the region of high protein accumulations transmembrane pressure differences between the lumen and ECS are virtually zero, creating a stagnant, motion-free region. The transient and steady-state radially-averaged axial protein concentration results from the one-dimensional and two-dimensional models are compared in Figure 5.6. The steadystate ECS profile was assumed to be reached at 1000 h as the maximum concentration change slope had decreased to less than 10-8 g/L-s. An ECS conductivity of 5x10-15 m2, a lumen flow rate of 500 mL/min and a low initial protein concentration of 5 g/L were selected for these computer model simulations in order to obtain relatively high radial ECS velocities. The results indicate that radial protein concentration gradients as well as radial ECS pressure gradients can be neglected from the analysis. Radial concentration variations at any axial location were not found for the 6 significant figures printed. To the same number of significant figures, the one-dimensional and two-dimensional models computed identical 55 radially-averaged axial concentration profile results. Hence the two-dimensional simulations, which required about two orders-of-magnitude more computation time, were unnecessary and the one-dimensional model was used to generate all further results. Figure 5.6 Transient axial protein concentration results from 2-D (-) and 1-D (.) model. The 2-D model results represent concentrations at R2 and R3, whereas the 1-D results represent radially-averaged concentrations. c30 = 10 g/L, Q 500 mL/min, k3= 5x10-15 M2, s = 1.0, L— 6X10-15 M. Insignificant radial protein concentration gradients in an empty ECS were also predicted by Taylor et al. (1993) and Patkar et al. (submitted). Patkar simplified Taylor's twodimensional model to a one-dimensional one and then compared the predictions of the two models. This comparison yielded a similar outcome, namely that the one-dimensional empty ECS model predicted very similar transient and steady-state ECS concentration profiles as the fully two-dimensional model. The radial Peclet numbers for an empty ECS HIFBR are so low (Pe << 1), that diffusive protein transport in the radial direction always dominates under 56 normal operating conditions. Axial and radial velocities in a packed ECS can be several orders-of-magnitude lower than in an empty ECS. Therefore, it is not surprising that the oneand two-dimensional models developed here for the latter case, yield such closely similar protein distribution results. 5.5.1 Comparison of Experimental Data with Model Results To validate the mathematical model developed here to describe transient protein transport in the ECS of HFBRs, a packed ECS was simulated experimentally by filling the ECS of gambro reactors with agarose/protein gels and then pumping a unidirectional flow of water through the fibre lumina. Several protein polarization experiments with different BSA loadings were performed. At defined times the reactors were frozen in liquid nitrogen and the axial ECS protein concentrations analyzed. The measured radially-averaged ECS protein concentration profiles were then compared with the model predictions. The protein polarization experiments are summarized in Table 5.4. In the first two experiments the fluid was not cleared from the lumen-side and therefore the ECS protein concentration calculations accounted for the dilution by the protein-free water in the lumen. In the other experiments the lumen fluid was first removed before freezing the reactor in liquid nitrogen. This procedure avoided cracks in the cartridge shell which complicated sawing of frozen reactors. In the experiments whose results are shown in Figures 5.8 and 5.11, a trace amount of azoalbumin (red-coloured BSA derivative) was added to the gel/BSA mixture to allow visual observation of the protein polarization in the HFBR. The red colour intensity increased at the downstream end of the reactor during continuous lumen flow, whereas gradual protein leakage from the upstream ECS manifold (Figure 5.1) resulted in red streaks along the cartridge wall. This latter material polarized less rapidly because protein transport by 57 convection and/or diffusion from the manifold region occurred very slowly. However, after the upstream manifold was completely cleared of protein, the red streaks disappeared and radial variations in the reactor were no longer apparent. The greatest red colour intensity was observed at the downstream end where the highest protein concentrations were measured. In the visually clear upstream regions very low protein concentrations were detected. Table 5.4 Summary of experimental protein polarization results. Protein Protein concentration concentration recovered added (g/1-) Protein Recovery Gel Time of concentration Recycle flow (w/w) (/0) (g/L) (/0) Lumen flow rate (h) (mL/min) 20 20 100 3 215.6 492 20 19.3 97 3 146.3 544 28.6 27.0 94 2 331.5 512 30 26.9 90 2 380.8 447 5 4.7 94 2 329.2 484 The sedimentation of protein solutions in the empty ECS of horizontally-operated HTBR cartridges was reported by Piret et al. (1991). In the experiments where azoalbumin was added to gel-filled HFERs (Figures 5.8 and 5.11), vertical concentration gradients due to solution sedimentation were not visibly noticeable. To evaluate the importance of the buoyant forces relative to the viscous resistance acting on the concentrated BSA solutions encountered near the downstream end of the HIFBR a Grashof number was calculated from: 58 e g p Ap Gr = ^ 2 (5 1) where dpore is the pore diameter of the gel, p the fluid density, g is the gravitational constant, Ap is the difference between the protein containing and protein free fluid densities and 1i is the viscosity of the solution. Using an average gel pore diameter of 1.0 i_tm (Bassi et al., 1987) and values listed in Table 3.1 for a solution containing 50 g/L of BSA, a Grashof number of <10-5 was calculated. Since this Grashof number was much less than unity, significant sedimentation of concentrated protein solutions was not expected. Computer modelling results using the measured agarose/protein (2%/2%) gel conductivity of 7.5x10-16 m2 did not agree with observed ECS protein profiles obtained for the first run. The computed ECS velocities were too low to cause a significant change from the initial uniform protein profile after 216 h (Figure 5.7). The gel structure in the conductivity measurement experiments was probably not the same as that in the cartridge. This may be explained by incomplete gel attachment at the fibre walls, which could have resulted in higher average ECS hydraulic conductivities. Since direct gel conductivity measurements on a gel filled HFBR cartridge were not possible, an apparent ECS conductivity of 5x10-15 m2 was obtained by forcing the transient computer model predictions to fit the experimental results of Figure 5.10. Comparisons for the other experimental runs (Figures 5.8 - 5.11) showed that, if a constant ECS conductivity of 5x10'5 m2was assumed, the computer model simulations adequately described all the measured transient ECS protein concentration profiles. The measured protein concentrations represent the average concentration of each cut section. In Figure 5.7, where the first section was about 10 cm long, it is obvious that the plotted concentration did not necessarily describe the true axial concentration profile. This experiment was repeated and more axial sections were analyzed to obtain the more detailed axial protein distribution shown in Figure 5.8, which revealed an improved agreement between the experimental and theoretical results. 59 .20 Axial Distance (cm) Figure 5.7 Experimentally obtained radially-averaged axial ECS protein concentrations versus corn puter modelling results. c30= 20.2 giL, Q = 492 mL/min, time = 215.6 h, 3% gel, k3= 5x10-15 and 7.5x10-16 m2, E = 1.0, Lp= 6x10-15m. Solid line: model;^experimental data, dashed line: section length. 50- 0 0^ 15 Axial Distance (cm) Figure 5.8 Experimentally obtained radially-averaged axial ECS protein concentrations versus computer modelling results. c30 = 19.3 g/L, Q = 544 mL/min, time = 116.3 h, 3% gel, k3= 5x10-15 m2, E = 1.0 Lo= 6x10-15m. Solid line: model;.: experimental data. 60 Higher initial protein loading experiments were performed (Figures 5.9 and 5.10) to test the influence of increased osmotic pressures on the protein redistribution behaviour. The protein transport was significantly slower and, after 14 - 16 days, nearly linear axial concentration profiles were measured. Again, the model provided a good description of the measured radially-averaged axial protein distributions in the ECS. To test the model at a lower initial protein concentration, an experiment was performed at 5 g/L BSA (Figure 5.11). As expected, the decreased osmotic effects increased the ECS velocities and protein was transported more rapidly to the downstream end of the reactor. The measured radiallyaveraged axial protein profile was no longer linear, but still agreed reasonably well with the model prediction. Computer simulation results were able to describe the axial protein concentration profiles measured in multi-fibre HFBRs. However, to account for an apparently altered gel structure in the reactors, it was necessary to assign a gel conductivity in the computer model that was about 10 times higher than the measured value. The experimental method of using an average protein concentration for a cut section, assumed that the protein was distributed uniformly in the radial direction. This assumption did not always hold since it was observed that higher azoalbumin concentrations near the cartridge wall existed in some regions, whereas the rest of the cross-section appeared to be clear of protein. By looking at the cut surface of the HFBR cross-sections, it was confirmed that the dye was mainly located near the cartridge walls. These higher protein concentrations could be due to the ECS manifold ring at the upstream end of the reactor, which released protein along its circumference into the outer regions of the fibre bundle. The proteins inside the multi-fibre bundle were more rapidly cleared from the upstream space by the secondary ECS convective flow. 61 50H T..;`) c 0 CO C 403020 - (i)^100 0 10 ^ 1'5 Axial Distance (cm) Figure 5.9 Experimentally obtained radially-averaged axial ECS protein concentrations versus computer modelling results. c30= 27.5 g/L, Q = 512 mL/min, time = 331.5 h, 2% gel, k3= 5x10-15m2, c = 1.0, Lp= 6X10-15m. Solid line: model;.: experimental data. 0 115 Axial Distance (cm) Figure 5.10 Experimentally obtained radially-averaged axial ECS protein concentrations versus computer modelling results. c30= 26.9 g/L, Q = 447 mL/min, time = 380.8 h, 2% gel, k3= 5x10-15m2, 6 = 1.0, Lp= 6x10'5 m. Solid line: model at 380.8 h; dotted line: steady-state result at 1888 h; •: experimental data. 62 Figure 5.11 Experimentally obtained radially-averaged axial ECS protein concentrations versus computer modelling results. c30= 4.7 g/L, Q = 484 mL/min, time = 329.2 h, 2% gel, k3= 5x10-15m2, E = 1.0, Lp= 6X10-15m. Solid line: model;.: experimental data. 5.5.2 Model Predictions The initial and steady-state ECS and lumen effective pressure distributions for a packed ECS (k3= 10'3m2) and for an empty ECS (k3= 10-9m2) are compared in Figures 5.12 and 5.13, respectively, for an initial BSA concentration of 10 g/L and a lumen flow rate of 500 mL/min. A packed bed ECS conductivity of 1043 m2 was chosen for Figure 5.12 since, for lower ECS conductivities the lumen and ECS pressures were no longer distinguishable. Due to the initially uniform protein distribution, the ECS pressure profiles for both conductivities were similar to the protein-free situation discussed in Section 5.4. At steady-state the proteins were polarized at the downstream end of the reactors, and the points where lumen and ECS pressures intersected were shifted upstream with decreasing ECS conductivities. 63 Axial Distance (cm) Figure 5.12 Lumen and ECS total pressures (P3 - 113) as function of axial distance at steady-state for a packed ECS conductivity. c30= 10 g/L, 0 = 500 mL/min, k3 = 10-131112, E = 1.0, Lp= 6X10-15 M. Figure 5.13 Lumen and ECS total pressures (P3 - 1-13) as function of axial distance at steady-state for an empty ECS. c30= 10 g/L, Q = 500 mL/min, k3 = 10-91112, L0= 6X10-15 M. 64 The transient protein concentrations and ECS velocity profiles for an empty ECS are plotted in Figure 5.14 and the associated initial and steady-state transmembrane velocities in Figure 5.15. Protein accumulation at the downstream end decreased transmembrane pressure differences and resulted in decreased ECS velocities. After steady-state was achieved, the secondary ECS flow was restricted to a region of about 75% of the total reactor length. Low BSA diffusivity resulted in a sharp concentration gradient between the downstream polarized and the upstream protein-free regions (Taylor et al., 1993; Patkar et al., submitted). The packed ECS conductivity of 10-13m2 required lower protein concentrations to counter-balance the lumen pressures since the transmembrane pressures were decreased. Hence, at steady-state, a larger ECS region contained proteins. For example, at steady-state, the proteins in the packed ECS were distributed over about 50% of the total fibre length (Figure 5.16), whereas in the empty ECS it was only about 25% (Figure 5.14). Also, due to the decreased axial ECS velocities (Figure 5.16), the time required to create the steady-state downstream polarized region was about 20 times longer. Diffusive protein transport, which did not influence the results in the empty ECS case, became more important in the packed ECS because of the greatly reduced ECS velocities. Hence the protein concentration gradient in the protein-containing regions was also determined by the protein diffusivity. The corresponding initial and steady-state radial transmembrane velocities are plotted in Figure 5.17. Whereas in the empty ECS case the maximum axial ECS (Figure 5.14) and transmembrane velocities (Figure 5.15) decreased significantly at steady-state, these maximum velocities changed little in the packed ECS case (Figure 5.16 and 5.17) A smoother, more spread out, steady-state transmembrane velocity profile, compared to Figure 5.15, separated the upstream protein depleted and downstream protein polarized regions (Figure 5.17), where low transmembrane fluxes still existed. Note, that the packed ECS conductivity used to generate Figures 5.16 and 5.17 was reduced back to its standard value of 5x1 5 2 65 Figure 5.14 Transient radially-averaged axial ECS protein distribution and axial velocities for an empty ECS. c30 = 10 g/L, Q = 500 mL/min, k3 = 109m2, L= 6x10-15 m. Solid line: concentration; dotted line: velocity. 30 iT) 40 20 s3) t ç30 0 8 10 steady-state 3 3 CT F13 CD 0^CD 20 -10 initial 10 Ea -30 0 0 5^10^15 20 Axial Distance (cm) Figure 5.15 Steady-state radially-averaged axial transient BSA concentrations and transmembrane velocities for an empty ECS, c30 = 10 g/L, 0 = 500 mL/min, k3 = 10-9 m2, Lp= 6x10-15 m. Solid line: concentration; dotted line: velocity. 66 35 300 30 -250 25 - -200 0 20 -150 es 5 15-1 100 4 -50 •^'0 10^15^20 Axial Distance (cm) Figure 5.16 Transient radially-averaged axial protein distribution and axial velocities for a packed ECS, c30 = 10 g/L, Q = 500 mL/min, k3 = 5x10'5 m2, c = 1.0, L— 6x10'5 m. Solid line: concentration; dotted line: velocity. Figure 5.17 Radially-averaged BSA concentrations and transmembrane velocities for a packed ECS, c30= 10 g/L, Q = 500 mL/min, k3= 5x10-15 m2, c = 1.0, L = 6x10'5 m. Solid line: concentration, dotted line: velocity. 67 The viscosity of BSA solutions can vary significantly from the viscosity of water as the BSA concentration is increased. In the model equations, the ECS fluid was treated as water and viscosity changes due to the presence of proteins were neglected. The computer program for the one-dimensional case was altered such that a concentration dependent viscosity was calculated to investigate the impact of viscosity changes due to the presence of BSA. Therefore Equation 3.32, governing the pressure distribution in the ECS, was rewritten as follows: d (1 d P3), ^2R1^L p d z 113^ dz (f? -1=0 k3 (pi p3 / 3 ) (5.2) where i_13 is the local viscosity of the ECS protein solution and 1..ti is the constant viscosity of the lumen and membrane fluid (water). Radially-averaged BSA concentrations were used to calculate the local ECS viscosity 43 using Equation 3.27a, developed by Tanford et al. (1956). After the axial protein distribution was recalculated at each new time-step, the associated ECS viscosities were updated. Since the viscous influence of the BSA solutions was expected to be more significant for higher protein loadings, an initial protein concentration of 20 g/L was selected for the model comparisons. The results were compared to those obtained using the earlier one-dimensional model, which did not include effects of variable viscosity. For both the empty and packed ECS cases, transient and steady-state radially-averaged axial protein profiles were investigated. The results obtained for the empty ECS case are shown in Figure 5.18. The concentration profiles are nearly identical to the cases where viscous effects were neglected at all times, indicating that the variations in protein solution viscosity had little effect on the behaviour of the system. 68 60 steady-state O h 1h 5^10^15 ^ 20 Axial Distance (cm) Figure 5.18 Transient radially-averaged ECS protein concentrations for an empty ECS, including viscosity variations. Solid line: constant viscosity, dotted line: variable viscosity. c30= 20 g/L, Q = 500 mL/min, k3= 109 m2, L = 6x10 15m. - steady state :LT 50 - CY) 500 h 40 Co -E. 30 U) (.) O h o 20 U) 10 CO o o 5^10^15 ^ 20 Axial Distance (cm) Figure 5.19 Transient radially-averaged ECS protein concentrations for a packed ECS, including viscosity variations. Solid line: constant viscosity, dotted line: variable viscosity. c30 = 20 g/L, Q = 500 mL/min, k3 = 5x10-15 m2, e = 1.0, Lp= 6x10-15m. 69 The radially-averaged protein concentration profiles for a packed ECS are shown in Figure 5.19. Because the flow resistance of the packed ECS plays a far more important role in controlling HIBR hydrodynamics than in the empty ECS case, the influence of variable viscosity was more pronounced. The transient protein concentration profiles lagged slightly behind the ones obtained when viscous effects were not considered. The steady-state criterion was satisfied after 1671 h in the latter case while an additional 95 h was required by the former to achieve the same steady-state concentration profile. In general, however, the results showed that ECS viscosity changes due to the presence of proteins did not significantly alter the transient concentration profiles. Also, at steadystate, the constant and variable viscosity models predicted nearly identical concentration profiles for both the empty and packed ECS cases. 280 0 .773, 60 4E- a) ° 40 0 0 ^ 0 5^10^15 ^ 20 Axial Distance (cm) Figure 5.20 Steady-state radially-averaged BSA concentrations for c30 = 5, 15, 30, 50 and 75 g/L, 0 = 500 mL/min, L= ^Solid line: packed ECS, k3 = 5x10-15 m2, c = 1.0; dotted line: empty ECS, k3 = 10-9 m2. 70 Radially-averaged steady-state ECS protein profiles for a packed ECS reactor having a conductivity of 5x10'5 m2and different initial protein concentrations are plotted in Figure 5.20. With increasing initial protein concentration, the protein profile extended over larger axial ECS regions. For protein loadings higher than 30 g/L the proteins were distributed throughout the entire ECS. The empty ECS radially-averaged protein concentration profiles obtained under otherwise identical conditions are also shown in Figure 5.20. For the empty ECS case, a minimal protein concentration of about 30 giL was also required to achieve complete protein distribution in the ECS. This result was somewhat surprising, since it was expected that a higher initial ECS protein loading would be required to achieve this result for the empty ECS case. Also, the steady-state profiles for higher initial concentrations were very similar to the ones obtained from the packed ECS model. At high initial protein concentrations (>30 g/L), sufficient protein is loaded in the ECS for osmotic pressures to almost exactly counter-balance the lumen pressure variation, which is the same in both the empty and packed ECS situations. The time required to reach steady-state as a function of initial protein concentration is plotted in Figure 5.21. With increasing protein loading the time required to reach steady-state initially increased, due to the reduced ECS velocities, caused by the osmotic back-pressure. However, above a critical protein concentration of about 25 g/L, the time began to decrease again, presumably because the proteins did not have to travel as far to reach their final steadystate distribution. For the empty ECS case the times are significantly lower, since axial ECS velocities are higher due to the higher ECS conductivity. Closer observation of Figure 5.20 reveals that, for the protein profiles obtained for initial concentrations below 30 g/L, true steady-state was not reached yet. Computer modelling results showed that about twice the time was required to avoid the low concentration tails in front of the downstream polarized protein profiles. By lowering the steady-state criterion these tails disappeared. But since these new protein concentration profiles were otherwise 71 very similar, the standard steady-state criterion^10-8 g/L-s) was maintained to generate all the results shown in Figure 5.21. Figure 5.21 Time to reach steady-state as a function of initial BSA concentrations. Q = 500 mL/min, Ic3 = 5x10-15 m2, c = 1.0, Lp = 6x10-15m. Solid line: packed ECS, dotted line: empty ECS The effect of ECS conductivity on the steady-state protein profiles is shown in Figure 5.22. For an empty ECS conductivity of 109m2, a highly polarized protein region existed at the downstream end of the reactor. With decreasing ECS conductivity the transmembrane pressure differences decreased and hence lower protein concentrations were required to shut off the transmembrane flows. Also, at reduced axial ECS velocities, diffusive transport became more dominant and the concentration gradients became more gradual. With decreasing ECS conductivities the time required to reach steady-state increased and, due to 72 low axial ECS velocities, a nearly linear concentration profile was obtained at an ECS conductivity of 10'7m2. Figure 5.22 Steady-state radially-averaged protein distributions at ECS conductivities of k3 = 10-9, 5x10-15, 10-16 and 10-17m2, e = 1.0, C30 10 g/L, Q = 500 mUmin, Lo= 6x10 15m. - The maximum average axial ECS Peclet number exists at the beginning, when proteins are uniformly distributed in the ECS. Based on the length of the fibre L, diffusivity of the protein D and the velocity obtained by averaging the ECS axial velocity over the entire ECS volume, this Peclet number was calculated as follows: Pe = u L ^(5.2) 3^ This maximum average axial ECS Peclet number, for different ECS conductivities is shown in Table 5.5. In the empty ECS case (1c3 10-9 m2) the high Peclet number shows that convective transport dominates. As the ECS conductivities decrease, the axial ECS velocities ^ 73 decreased and hence the Peclet numbers are lower. For the lowest ECS conductivity of 10-17 m2, convective transport no longer dominates. Table 5.5 Maximum average axial ECS Peclet numbers. ECS conductivity^Peclet number k3^ Pe [m2]^ [-] 10-9^16200 5x10-15^216 1047^ 0.4 In Figure 5.23 the transient radially-averaged axial protein concentration and axial ECS velocities are plotted for the tissue-like packed ECS conductivity of 10-17 m2. At steady-state the axial ECS velocities do not significantly differ from their initial values. The steady-state profiles were obtained after about 10 years at which time only a slight protein gradient existed. Normally HFBRs operate for about 3 months. Over such times, at low ECS conductivities, protein polarization should therefore be negligible. The BSA diffusivity used in the computer simulations was considered to be concentration independent (van den Berg and Smolders, 1987). The modelling results were obtained at protein and fluid properties evaluated at room temperature (20 °C). The increased HFBR operating temperatures used for mammalian cell culture (37 °C) will increase the diffusivities of proteins. As a result, axial diffusive transport will increase further and protein polarization may be even less significant. 74 Axial Distance (cm) Figure 5.23 Transient radially-averaged BSA concentrations and ECS velocities for a tissue-like conductivity, c30= 10 g/L, Q= 500 mL/min, k3= 107m2, £ = 1.0, L= 6x10-15m. Solid line: concentration, dotted line: velocity. In the previous modelling the ECS porosity was taken to be 1. The influences of reduced ECS volumes due to the presence of impermeable cells in the ECS were investigated. The results for a low ECS conductivity of 10-17m2 are shown in Figure 5.24. Decreased bed porosities increased the axial ECS velocities and this resulted in increased convective protein transport. Also, reduced bed porosities increased the downstream protein polarization. The shortest period of time to reach steady-state was required for the lowest ECS porosity. 75 Axial Distance (cm) Figure 5.24 Steady-state radially-averaged protein distributions in fluid phase at ECS porosities of -= 1.0, 0.5, 0. 2 and 0.1, k3 = 10-17m2, c30 = 10 g/L, Q = 500 mL/min, Lp= 6X10-15 M. 76 CHAPTER 6 Conclusions In this work, the transport of osmotically-active, high molecular weight proteins entrapped in the ECS of BFBRs was studied. One- and two-dimensional mathematical models were developed to describe the motion of fluid and protein in HFBRs whose ECS was packed by cells or gel. Numerical methods were required to solve the coupled differential equations governing the transient ECS protein distribution and quasi-steady hydrodynamics in the lumen and ECS. The packed bed of cells was simulated experimentally by filling the ECS with liquid agarose solution containing protein. Upon cooling, a porous gel was formed having a sufficiently low hydraulic conductivity that Darcy's law was applicable. Therefore, monthlong experiments, required to grow mammalian cell packed HFBRs, were not needed. The computer modelling results indicated that radial ECS pressure gradients were negligible and hence, the simplified one-dimensional model was sufficient for all HFBR operating conditions of interest. The latter computer program required about two orders-ofmagnitude less time to produce the same results as the two-dimensional model. 6.1 Experimental Versus Model Results Axial ECS velocities were mainly governed by the hydraulic conductivity of the gel. When the gel conductivity obtained by the falling head method (Bear, 1972) was used, the model was unable to describe the experimental axial protein profiles measures in the ECS of gel filled HFBRs. It was believed that the gel did not completely fill the complex ECS volume. Thus, a higher conductivity was obtained by fitting the model predictions to one set of experimental 77 data. The modelling results indicated that, with an effective ECS conductivity of 5x1015 m2, all of the experimentally obtained axial protein concentrations profiles could be adequately described. The primary differences between model and experimental results occurred near the ends of the reactor and were attributed to the circumferential ECS manifolds, which exchange fluid and protein only very slowly with the fibre bundle region of the ECS. 6.2 Implications for Hollow Fibre Bioreactor Operation Mammalian cells entrapped in the ECS of a HFBR can eventually fill this space to produce very highly packed conditions with low hydraulic conductivities. The hydraulic conductivity of hybridoma cell lines can be roughly approximated by those for packed beds of red blood cells; values of which have been reported by Zydney et al. (1986) and by packed beds of spheres whose conductivities can be described by the Kozeny-Carman equation (Zydney et al., 1986). Surface dependent cell lines like BHK (baby hamster kidney) form an extracellular matrix around the cells and conductivities are expected to be like those of tissues, which are orders-of-magnitude lower. Modelling results indicated that ECS velocities, responsible for the convective protein transport, were so low at packed ECS conductivities that, practically no protein polarization occurred on the time scale of HFBR operation. The hydrodynamic resistance in the ECS became so dominant that the membrane hydraulic permeability no longer affected the ECS flow. The ECS pressure distribution, nearly constant in an empty ECS, was nearly identical to the lumen pressures. Hence transmembrane pressure differences, responsible for the leakage, flow decreased dramatically. Convective protein transport and associated osmotic effects influenced the ECS hydrodynamics much less since the protein profile changed little even after several months. 78 High BSA loadings at HFBR start-up can reduce downstream polarization of proteins (Taylor et al., 1993; Patkar et al, submitted). A uniform distribution of essential cell growth proteins (e.g., insulin) is desirable, especially upon cell inoculation, to ensure homogeneous cell growth. Modelling results showed that, after a high cell density is obtained, these high protein loadings are no longer required since convection no longer dominates the protein transport. Large product molecules (e.g., antibodies) entrapped in the ECS are harvested by opening the shell ports. The modelling results indicated that, in highly packed tissue-like cell beds, the ECS flow resistance is so high that practically no ECS leakage flow occurred. Consequently, after the ECS ports are opened it could take up to several weeks to recover the ECS fluid containing the desired product. Alternatively one could apply an external flow through one of the ECS ports to increase the flow of product-containing fluid from the other ECS port. However, to increase the harvesting flow significantly, very large pressures at the ECS side would be required. These high pressures could conceivably cause cell bed compression and hence further decreases in bed conductivities. Such decreases in packed cell bed conductivities with increasing pressure were reported by Zydney et al. (1986) for red blood cell filtration. The low ECS conductivities greatly reduced ECS fluid flow during product harvesting. To allow recovery of the ECS proteins the membrane size could be increased (i.e., microporous hollow fibre membranes), allowing the product to penetrate into the lumina. Alternatively, a dual circuit HFBR with two different membrane fibre cut-offs could be used to separate product harvesting from medium perfusion, but this would require a more complex reactor design. Local depletion of growth factors in highly metabolically active cell beds could cause cell starvation. Due to their high molecular weights these molecules cannot penetrate through the fibre membranes into the ECS and hence must be supplied from the ECS ports. Decreased ECS porosities increased the convective protein transport in the ECS even for a very low ECS 79 conductivity. As long as the ECS conductivity is not too low, convection can facilitate the transport into the metabolically active regions. But in highly packed ECS regions, convection no longer dominates. Transport of growth factors into these regions will now be governed by slow convective and diffusive transport; hence, cell growth may decrease or cell death may occur. 80 CHAPTER 7 Future Work 1) The fluid dynamics of HFBRs having hydraulic conductivities corresponding to packed beds of mammalian cells has not been thoroughly investigated. Ex situ techniques, such as the falling head method, did not provide measured conductivities representative of those found in HFBRs whose ECS is packed with the same porous medium. Due to the extracellular matrix formed around mammalian cells, conductivity calculations based on a spherical cell packing cannot be used to obtain the required estimates, especially in the geometrically complex ECS. A method should be designed to allow for accurate conductivity measurements in situ. 2) An on-line non-invasive flow measuring technique such as nuclear magnetic resonance could give experimental insight about the effects of osmotically-active proteins on ECS hydrodynamics. If this technique is sensitive enough, whole-reactor radial variations could also be investigated. 3)^The Krogh cylinder approximation assumes that the hydrodynamic and protein transport behaviour observed for a single fibre is representative of all fibres in the multi-fibre bundle which makes up the HFBR. A more sophisticated model could be developed to provide more realistic simulations of a multi-fibre HFBR for circumstances where the Krogh cylinder approximation is no longer applicable. Such a model could accommodate the existence of circumferential manifolds near the ends of the HFBR and would permit more realistic modelling of open-shell operations such as exist, for example, during product harvesting. This model could also handle the sedimentation of protein solutions that occur due to concentration-dependent density gradients. 81 4) During cell growth cell densities increase and tissue-like regions are formed that change the HFBR hydrodynamics. The current models could be re-written to include such time dependent local conductivity changes. 5) Smaller molecules such as insulin, which are essential for cell growth, can penetrate through the membrane. Protein polarization profiles need to be studied, both theoretically and experimentally, to understand the influence of membrane leakage on the behaviour of the system. For such studies, the impact of the simultaneous presence of high molecular weight proteins, such as BSA, could be included. 82 Nomenclature [m3 kg-1] [m6 kg-2] A3^3rd virial coefficient^ A cap^cross-sectional area of the capillary ^ [m2] [m2] A gel^cross-sectional area of the gel^ A mem^total area of the membrane in the HFBR^[m2] A 2^2nd virial coefficient^ c^concentration^ [kg m-3] c 30^uniform initial start concentration in ECS^[kg m-3] Z.-^radially-averaged axial concentration^[kg m-3] [m2 s-1] diffusivity^ D D.^diffusivity at a infinite diluted solution^[m2 s-1] dpore^pore diameter^ [m] 6^bed porosity^ F^body force vector^ gravitational constant^ g H [N] [m s-2] h^height of water in the capillary^ [m] h0^height of constant water head^ [m] I.^modified Bessel function, first kind, order zero ^H II^modified Bessel function, first kind, order one^H membrane flux^ [m3 s-1 m2] J [m2] k^conductivity^ pressure drop parameter (Park and Chang, 1986)^H K K constant for protein viscosity Eqn. 3.27a ^H IC^convective hindrance factor^ H Kd^diffusivity hindrance factor ^ H K.^modified Bessel function, second kind, order zero^H K1^modified Bessel function, second kind, order one^H L reactor length^ Lp^membrane permeability^ M^molecular weight^ nis^molar salt concentration^ N P O number of fibres in the reactor^ pressure^ [m] [m] [kg moles-1] [g L-1] H [Pa] volumetric flow rate through the HFBR ^[m3 s-1] 83 Qi^volumetric flow rate through a single fibre ^[m3 s-11J [m3 s-1] Q3^volumetric flow rate^ aleak^ average volumetric leakage flow rate from the lumen ^[m3 s-1] r^radial distance^ R^gas constant^ R1^lumen radius^ R2^outer membrane radius^ R3^Krogh radius^ Rc^cartridge radius^ S particle surface area^ I^time^ T^absolute temperature^ u^axial velocity^ Ft^radially-averaged axial velocity^ , [m] [J moles-1 K-1] [mi [m] [m] kill [m2] [s] [K] [m s-1] [m s-1] u3^ axial ECS velocity averaged over the entire ECS ^[m s-1] v^radial velocity^ [m s-1] ^ ^ [m s-1] [m3] Z^protein charge number^ [ri] [-] velocity vector^ particle volume^ z^axial distance^ Abbreviations BSA^bovine serum albumin ECS^extracapillary space EIFBR^hollow fibre bioreactor ODE^ordinary differential equation PDE^partial differential equation 84 Subscripts and Superscripts 1^lumen region 2^membrane region 3^ECS region averaged quantity i^index for axial grid location j^index for radial grid location in^fibre inlet out^fibre outlet initial situation o Greek Letters 5^thickness of the gel^ porosity^ £ 1-1^viscosity^ II.^viscosity at infinite diluted solution^ [Il]^intrinsic viscosity^ it^3.14159....^ [In] [-] [Pa-s] [Pa-s] [L g-1] p^density^ [-] [Pa] [kg m-3] a^protein reflection coefficient^ kinematic viscosity^ u [-] [m2 s-1] T^flow potential^ AP^pressure difference between the lumen [Pa] H^osmotic pressure^ inlet and ECS outlet^ Dimensionless Numbers Gr^Grashof number Pe^Peclet number Re^Reynold's number [Pa] 85 References Anderson, J.L., Rauh, F., Morales, A., "Particle Diffusion as a Function of Concentration and Ionic Strength", J. Phys. Chem., 82, 608-6016 (1978). Anderson, D.A., Tannehill, J.C., Plechter, R.H., Computational Fluid Mechanics and Heat Transfer, Hemisphere, New York, 1984. Apelblat, A., Katzir-Katchalsky, A., Siberberg, A., "A Mathematical Analysis of Capillarytissue Fluid Exchange", Biorehol , 11, 1-49 (1974). Atkinson, B., Mavituna, F., Biochemical Engineering and Biotechnology Handbook, 2nd Edition, Nature Press, 1991. Bassi, A.S., Rohani, S., Macdonald, D.G., "Measurement of Effective Diffusivities of Lactose and Lactic Acid in 3% Agarose Gel Membrane", Biotechnol. Bioeng., 30, 794-797 (1987). Bear, J., Dynamics of Fluids in Porous Media, Dover Publications, New York, 1972. Belfort, G., "Membranes and Bioreactors: A Technical Challenge in Biotechnology", Biotechnol. Bioeng., 33, 1047-1065 (1989). Berman, A.S., "Laminar Flow Channels with Porous Walls", J. Appl. Phys., 24, 1232-1235 (1953). Broek, A.P., Bargeman, D., Sprengers, E.D., Smolders, C.A., "Characterization of Hemophan Hemodialysis Membranes by Thermoporometry", Intl. J. Artif. Org., 15, 25-128 (1992). Bruining, W.J., "A General Description of Flows and Pressures in HF-Membrane Modules", Chem. Eng. Sci., 44, 1441-1447 (1989). Carman, P.C., "Fluid Flow trough a Granular Bed", Trans. Inst. Chem. Eng. London, 15, 150156 (1937). Chen, T.L., Humphery, A.E., "Estimation of Critical Particle Diameters For Optimal Respiration of Gel Entrapped and/or Pelletized Microbial Cells", Biotechnol. Letts., 10, 699 (1988). Collins, R.E., Flow of Fluids Through Porous Materials, Reinhold Publishing, New York, 1961. Cooney, D.O., Kim, S., Davis, J., "Analyses of Mass Transfer in Hemodialyzers for Laminar Blood Flow and Homogeneous Dialysate", Chem. Eng. Sc., 29,1731-1738 (1974). Darcy, H., Les Fontaines Publiques de la Ville de Dijon, Dalmont, Paris, 1856 86 Dean, P.D.G., Johnson, W.S, Middle, F.A., Affinity Chromatography, a Practical Approach, IRL Press, New York, 1985. Gin, H., Dupuy, B., Caix, J., Baquey, Ch., Ducassou, D., "In vitro Diffusion in Polyacrylamide Embedded Agarose Microbeads", J. Microencaps., 7, 12-23 (1990). Hammer, B.E., Heath, C. A. , Mirer, S.D., Belfort, G., "Quantitative Flow Measurements in Bioreactors by Nuclear Magnetic Resonance Imaging", BiotechnoL Bioeng., 8, 327 (1990). Heath, C., Belfort, G., "Immobilization of Suspended Mammalian Cells: Analysis of Hollow Fiber and Microcapsule Bioreactors", Adv. in Biochem. Eng./Biotechnol, 34, 1-31 (1987). Heath, C., Belfort, G., Hammer, B.E., Mirer, S.D., Pimbley, J.M., "Magnetic Resonance Imaging and Modeling of Flow in Hollow-Fibre Bioreactors", A.LCh.E. 1, 4, 547-558 (1990). Hjerten, S., "Molecular Sieve Chromatography on Polyacrylamide Gels, Prepared According to a Simplified Method", Arch. Biochem. Biophys., 1, 147-151 (1962). Humphrey, A., Aiba, S., Millis, N., Biochmechanical Engineering, Univ. Tokyo Press, Tokyo, Japan, 1985. Kelsey, L.J., Pillarella, M. R., Zydney, A. L., "Theoretical Analysis of Convective Flow Profiles in a Hollow-Fiber Membrane Bioreactor", Chem. Eng. Sci., 11, 211 (1990). Klein, E., Autian, J., Bower, J.D., Buffaloe, G., Centella, L.J., Colton, C.K., Darby, T.D., Farrell, P.C., Holland, F.F., Kennedy, R.S., Lipps, B., Mason, R., Karl, D.N., Villarroel, F., Wathen, R.L., Evaluation of Hemodialyzer and Dialysis Membranes, Report of a study group for the artificial kidney - chronic urea program, U.S. Department of Health Service, DHEW Publication, Washington D.C., 1977. Kleinstreuer, C., Poweigha, T., "Modeling and Simulation of Bioreactor Process Dynamics", Adv. in Biochem. Eng./Biotechnot, 30, 91 (1984). Kleinstreuer, C., Agarwal, S.S., "Analysis and Simulation of Hollow-Fiber Bioreactor Dynamics", Biotechnot Bioeng., 28, 1233-1240 (1986). Knatzek, R.A., Guillino, P. M., Kohler, P. 0., Dedrick, R. L., "Cell Culture on Capillaries: an Approach to Tissue Growth in vitro", Science, 178, 65-67 (1972). Knatzek, R.A., Kohler, P. 0., Guillino, P. M.," Hormone Production by Cells Grown in vitro on Artificial Capillaries", Exptl. Cell Res., 84, 251-254 (1974). Kozeny, J., Hydraulic, Springer, Vienna, 1953. 87 Krogh, A., "The Number and Distribution of Capillaries in Muscles with Calculations of the Oxygen Pressure Head Necessary to Supplying the Tissue", J. Physiot , 52, 79 (1919) Levick, J.R., "Flow Through Interstitium and other Fibrous Matrices", Quart. J. Experim. PhysioL, 72, 409-437 (1987). Libicki, S.B., Salmon, P.M., Robertson, C.R., "The Effect Diffuse Permeability of the Nonreacting Solute in Microbial Cell Aggregates", Biotechnol. Bioeng., 32, 64 (1986). Libicki, S.B., "Transport of Gases and Liquids through Dense Microbial Cell Aggregates Cultured within Hollow Fibre Membrane Bioreactors", Ph.D. Thesis, Stanford (1986) Mulder, M., Basic Principles of Membrane Technology, Kluwer Academic Publisher, Dordrecht, 1991. Pallmann, H., Deuel, H., "Ober die Wasserdurchlassigkeit von Hydrogelen", Experimenta, 1, 325-326 (1945). Park, J.K., Chang, H.N., "Flow Distribution in the Fiber Lumen Side of a Hollow-Fiber Module", A.I.Ch.E. J, 12, 1937-1947 (1986). Patankar, S.V., Numerical Heat Transfer and Fluid Flow, Hemisphere Publ. Co., New York, 1980. Patkar, A.Y., Bowen, B.D., Piret, J.M., "Protein Adsorption to Polysulfone Membranes in Hollow Fibre Bioreactors used for Serum-free Mammalian Cell Culture", BiotechnoL Bioeng., in press, (1993). Patkar, A.Y., Koska, J., Taylor, D.G., Bowen, B.D., Piret, J.M., "Experimental and Modeling Study of Hollow Fiber Bioreactor Protein Distribution", submitted. Peaceman, D.W., Rachford, H.H., "The Numerical Solution of Parabolic and Elliptic Differential Equations", J. Soc. Ind. AppL Math., 3, 28-41 (1955). Pillarella, M.R., Zydney, A.L., "Theoretical Analysis of the Effect of Convective Flow on Solute Transport and Insulin Release in Hollow Fiber Bioartificial Pancreas", J. Biomech. Eng., 112, 220-228 (1990). Piret, J.M., Cooney, C.L., "Immobilized Mammalian Cell Cultivation in Hollow Fiber Bioreactors", Biotechnot Adv., 8, 762-783 (1990a). Piret, J.M., Cooney, C.L., "Mammalian Cell and Protein Distributions in Ultrafiltration Hollow Fiber Bioreactors", Biotechnot Bioeng., 36, 902-910 (1990b). Piret, J.M., Cooney, C.L., "Model of Oxygen Transport Limitations in HFBRs", Biotechnot Bioeng., 37, 80-92 (1991). 88 Piret, J.M., Devens, D.A., Cooney, CL., "Nutrient and Metabolite Gradients in Mammalian Cell Hollow Fiber Bioreactor", Can. J. Chem. Eng., 69, 421-426 (1991). Rony, P.R., "Multiphase Catalysis. II: Hollow-Fiber Catalysis", BiotechnoL Bioeng., 13, 431 (1971). Saffman, P.G., "On the Boundary Condition at the Surface of a Porous Medium", Stud. Appl. Math. , 50, 93-101 (1971). Salmon, P.M., Libicki, S.B., Robertson, CR., "A Theoretical Investigation of Convective Transport in the Hollow-Fiber Reactor", Chem. Eng. Comm., 66, 221-248 (1988). Sangani, A., "Transport Processes in Periodic Arrays of Spheres", Ph.D. Thesis, Stanford (1982). Sardonini, C.A., DiBiasio, D., "An Investigation of the Diffusion-Limited Growth of Animal Cells Around Single Hollow Fibers", Biotechnol Bioeng., 40, 1233-1242 (1988). Scatchard, G., Batchelder, A.C., Brown, A., "Preparation and Properties of Serum Plasma Proteins. VI. Osmotic Equilibria in Solutions of Serum Albumin and Chloride", J. Amer. Chem. Soc., 68, 2320-2329 (1946). Scheidegger, A.E., The Physics of Flow Through Porous Media, 2nd edition, University of Toronto Press, Toronto, 1960. Starling, E.H., "On the Absorption of Fluids from the Convective Tissue Space", 1 PhysioL , 19, 312-326 (1896). Swabb, E., Wei, J., Gullion, P. M., "Diffusion and Convection in Normal Neoplastic Tissues", Cancer Res., 34, 2814-2822 (1974). Tanford, C., Buzzell, J.G., "The Viscosity of Aqueous Solutions of Bovine Serum Albumin Between pH 4.3 and 10.5", 1 Phys. Chem., 60, 225 (1956). Tanford, C., Physical Chemistry of Macromolecules, John Wiley & Sons, New York 1961. Taylor, D.G., Bowen, B.D., Piret, JM., "Protein Polarization in Isotropic Membrane Hollow Fiber Bioreactors", A.I.Ch.E. J, in press (1993). Taylor, G.I., "A Model for the Boundary Condition of a Porous Medium", 1 Fluid Mech., 49, 319-326 (1971). Tharakan, J.P., Chau, P.C., "Operation and Pressure Distribution of Immobilized Cell Hollow Fiber Bioreactors", Biotechnot Bioeng., 28,1064-1071 (1986). Tharakan, J.P., Chau, P.C., "Modeling and Analysis of Radial Flow Mammalian Cell Culture", Biotechnot Bioeng., 29, 657-671 (1987). 89 Tombs, M.P., Peacocke, A.R. The Osmotic Pressure of Biological Macromolecules, Clarendon Press, Oxford, 1974. van den Berg, G.B., Smolders, G.B., "The Boundary-Layer Resistance for Unstirred Ultrafiltration. New Approach", J. Membrane Sci., 40, 149-172 (1989). van't Hoff, J.H., "Die Rolle des osmotischen Druckes in der Analogie zwischen LOsungen and Gasen", Z. Phys. Chem., 1, 481-508 (1887). Vilker, V.L., Colton, C.K., Smith, K.A., "The Osmotic Pressure of Concentrated Protein Solutions: Effect of Concentration and pH of Saline Solutions of Bovine Serum Albumin", Coll. Interface Sc., 79, 548-565 (1981a). Vilker, V.L., Colton, C.K., Smith, K.A., "Part II: Theoretical and Experimental Study of Albumin Ultrafiltered in an Unstirred Cell", J. Coll. Interface Sci., 27, 637-645 (1981b). Walker, J.M., Gingold, E.B., Molecular Biotechnology & Biotechnology, Royal Society of Chemistry, London, 2nd Edition, 1988. Waterland, L.R., Michaels, A.S., Robertson, C.R. ,"A Theoretical Model for Enzymatic Catalysis Using Asymmetric Hollow Fiber Membranes", A.I.Ch.E. 1. , 20, 50-59 (1974). Waterland, L.R., Michaels, A.S., Robertson, C.R., "Enzymatic Catalysis using Asymmetric Hollow Fiber Membranes", Chem. Eng. Comm., 2, 37-47 (1975). Webster, TA., Shuler, M.L. "Mathematical Models of Hollow-Fiber Enzyme Reactors", Biotechnol. Bioeng., 20, 1541-1556 (1978). Webster, IA., Shuler, M.L., "Whole-Cell Hollow Fiber Enzyme Reactor: Effectiveness Factors", Biotechnol. Bioeng., 21, 1725 (1979). Webster, I.A., Shuler, M.L., "The Whole Cell Hollow Fiber Reactor-III", Chem. Eng. Sci., 34, 1273 (1979). Webster, I. A., Shuler, M. L., "Whole-Cell Hollow-Fiber: Transient Substrate Concentration Profiles", Biotechnol. Bioeng., 22, 447-450 (1981). Wei, J., Russ, M.B., "Convection and Diffusion in Tissue and Tissue Cultures", J. Theor. Biol., 66, 775-787 (1981). Zydney, A.L., Saltzman, W. M., Clark, K.C., "Hydraulic Resistance of Red Cell Beds in an Unstirred Filtration Cell", Chem. Engng. Sci, 44, 147-159 (1986). 90 Appendix 1 Two-Dimensional Analytical Solution for HFBR Hydrodynamics in a Packed ECS, Protein Free System This is a typographically corrected version of the Apelblat et al. (1974) solution. In the absence of osmotically-active molecules, the pressure distributions in the fibre lumen, in the fibre membrane and in the ECS region can be calculated as follows: For the fibre lumen: pi(z) V= u= An^ D +G(z , ^ coskanz) n=1 an 1 dP, (1,2 R12) dz 1 d2 P ^ ^ (2r R,2 –t-3) 16p d z2 For the membrane and cell region ( i = 2, 3): L p(z,r) D+G(—+II3nR (anr) cos(anz)) 2 91 u tz,r ap, [I. az _ k k v,(z,r )= a P: IA a r Volumetric flow rate into the lumen: = icR 4 1 G Average volumetric leakage flow rate from the lumen: 4^co n R1 (117C) sin 8p,^n,1^n^2 QI,leak ^ GE A = pri n 7t g =^ An 13i, r4,1 an^" A 13 D=I t^n1 -FGZ an" 16k \ A =^SnICI VtnRi + WnI, (ar,R1 )) n^R1-3 92 2[(_on _11 13n a„L(A.„ +a„ 13) 13„ = W„ I JanR,)+S„Ko(a„Ri) E„ = K,(ar,R3) IjanR2) Kja„R, ManR3) D Io(a R2) ^K,(a„R.3)^ +K0(a„R2) Mat,R3) P—P G—^ out L—E R2 (anr) = W„Ijanr) + S.K.(anr) Ij a ^fl R^= K(anr) +K j^ o: nR3) ( S WT, /) a ^a„R2(D„II (a n R2) — Enlo (an R2 ) k k 2) 93 Appendix 2 One-Dimensional Analytical Solution for Hydrodynamics in a Packed ECS, Protein Free System In the absence of osmotically-active molecules the one-dimensional pressure distributions in the fibre lumen and ECS region can be calculated as follows: For the fibre lumen: Pi(z) = B, sinh(X4 +B2 cosh(kz) +B3 z +134 For the ECS: P3(z) =7 B, sinh(kz)+7 B2 cosh(2z)+B3 z +B4 Volumetric flow rate into the lumen: Tc R4 , = ^ 03, X, +B3) 81A Average volumetric leakage flow rate from the lumen: 4 L L (2, jeak = n R1 [A,B, (cosh(X— )– 1)+ B2 sinh (X — 8 t^ [ 2^2 Y= 1^R14 8k3 (R; – R22) ^ 94 , ^2R^L^16L p (R32 - R22) k 3^R3 sinh(A, L) B^ 1=^ cosh(X, L) — 1 2 B2 = Pow — P,,, sinh(X L) cosh(XL) ^ [sinh(XL)— X y Li — 1 cosh(A,L)— 1 B3 = — X y B, 95 Appendix 3 Experimental HFBR Protein Results Figure 5.7 Axial location of cut HFBR section [cm] 0.0^-^10.0 10.0^-^15.0 15.0^-^18.0 18.0^-^20.0 20.0^-^21.5 section weight water added (without HFBR shell) to section average measured BSA concentration in gel/water mixture [g] [mL] [g/L] 80.6 35.8 20.1 12.7 25.0 121.6 59. 53.9 39.7 57.0 0.88 4.38 4.48 4.79 7.06 25.0 34.8 23.3 14.7 15.5 13.5 14.2 14.0 20.5 30.8 42.2 52.0 30.3 30.8 30.0 29.9 45.1 35.5 0.66 2.62 1.69 1.80 3.13 3.80 4.10 3.95 8.70 21.3 10.6 11.7 11.4 10.3 11.3 11.2 12.9 12.1 8.8 24.5 90.6 59.8 60.9 61.2 59.5 62.2 54.1 60.4 75.4 46.8 97.3 0.12 0.96 1.48 2.12 2.54 3.65 4.32 4.65 4.00 4.85 8.95 Figure 5.8 0.0^-^2.0 2.0^-^7.0 7.0^-^10.0 10.0^-^12.0 12.0^-^14.0 14.0^-^16.0 16.0^18.0 18.0^20.0 20.0^21.5 Figure 5.9 0.0^2.2 2.2^-^3.9 3.9^6.0 6.0^8.1 8.1^9.9 9.9^12.0 12.0^13.7 13.7^15.8 15.8^17.6 17.6^18.9 18.9^-^21.5 96 Figure 5.10 Axial location of cut FIFER section section weight water added (without HFBR shell) to section [cm] 0.0^-^2.5 2.5^-^5.0 5.0^7.5 7.5^-^10.0 10.0^-^12.5 12.5^-^15.0 15.0^-^17.0 17.0^-^19.4 19.4^-^21.5 average measured BSA concentration in gel/water mixture [g] [mL] [g11-] 23.9 14.5 15.6 15.4 14.1 15.6 11.7 14.9 23.0 200 200 200 200 200 200 200 400 400 0.53 0.60 0.97 1.22 1.35 1.80 1.57 1.10 2.13 24.3 23.2 22.9 18.7 12.2 11.5 9.80 22.8 130 230 204 233 241 243 228 448 0.00 0.08 0.11 0.07 0.08 0.20 0.31 0.69 Figure 5.11 0.0^-^2.5 2.5^-^6.6 6.6^-^10.4 10.4^-^13.4 13.4^-^15.4 15.4^-^17.4 17.4^-^19.2 19.2^-^21.5 Sample Calculation for Gel BSA Concentration (fibre lumen was cleared from water) e.g. , Figure 5.11 last sample: length of axial cut AL = 21.5 - 19.2^2.3 cm weight of fibres^= 1.47 cm • 2.3 g/cm = 3.8 g weight of gel^= 22.8 g - 3.8 g^19 g volume added to gel/fibre mixture ^448.4 mL measured average concentration of BSA in gel/water mixture^0.69 g/L concentration of BSA in gel = ( 448.4 mL + 19 mL) / 19 mL • 0.69 g/L = 17.0 g/L 97 Experimental Gel Conductivity Results Area of the gel cross-section^A gel = 12 57 [cm2] Area of the capillary^kap = 4.63^[mm2] . 3% agarose gel time water height height of water head t [s] h ho [cm] [cm] [cm] 4310 7332 16805 1 2.3 3.5 17.1 17.1 17.1 1.95 1.95 1.95 2.4 7.8 9.3 10.4 15.6 17.2 18.3 18.3 18.3 18.3 18.3 18.3 2.15 2.15 2.15 2.15 2.15 2.15 1.2 1.6 2 10.7 10.7 10.7 2.05 2.05 2.05 0.3 0.5 0.6 2.15 2.5 2.7 9.9 9.9 9.9 9.9 9.9 9.9 2.3 2.3 2.3 2.3 2.3 2.3 1 2.1 4 5 9 10.5 13.2 19.6 19.6 19.6 19.6 19.6 19.6 19.6 2.5 2.5 2.5 2.5 2.5 2.5 2.5 axial gel thickness °pet 2% agarose gel 2820 12420 16380 19500 32700 78000 2% agarose gel 4620 6300 7860 2% agarose gel + 2% BSA 275 520 617 2700 3260 3615 1% agarose gel 111 272 632 842 2350 3320 6600 98 Appendix 4 The following pages lists the source-codes for the two-dimensional and one-dimensional computer models written in the FORTRAN language. Each program requires a separate input file in which all of the system parameters are specified. The over-relaxation parameters, required for the two-dimensional model are grid size and operating condition dependent; therefore, the parameters listed in the example input file are not necessarily the optima for different conditions. Also the COMMON block size can be easily changed (in both programs) to allow for larger axial and/or radial grids. Also, to avoid data loses due to run-time errors, individual files are opened for each transient solution (starting with the extension *.000 for the initial start-up data to a maximal *.099). In the two-dimensional model, two files are opened, one containing the axial and radial solutions, while the other contains only radially-averaged values. After the programs are initialized, an information file is written specifying all model parameters. In the one-dimensional model the protein distribution in an empty ECS can be generated by setting the ECS conductivity to 0. The program automatically calculates an ECS conductivity which is defined by geometric factors for the one-dimensional empty ECS situation. 99 Source Code for the Two-Dimensional Computer Model PROGRAM HFBR_Protein_2D * This program models the 2 dimensional transient protein concentrations in a single fibre HFBR. VERSION 4.1 last change : 6.6.93 by liirgen Koska Important variables used in the program • • • * • • • • • • • • • • • * • * • • • • • • • • • • • • * • • • • • • C_ECS(z,r)^2 dimensional field with ECS protein concentrations P_ECS(z,r)^2 dimensional field with ECS pressures P_LU (z)^I dimensional field with lumen pressures P_OS (z)^I dimensional field with osmotic ECS pressures X_^1 dimensional filed for axial grid location D_X^I dimensional field for axial grid spacing R_^I dimensional field for radial grid location D_R^I dimensional field for axial grid spacing RI^inner membrane radius R2^outer membrane radius R3^Krogh radius P_IN^lumen inlet pressure POUT^lumen outlet pressure PI^3.1415... P_EPS^accuracy for pressures in ECS and lumen NITER_P^max. 6 of iterations to reach PEPS ALPHAX^axial ECS pressure over-relaxation parameter ALPHAR^radial ECS pressure over-relaxation parameter DT^time step T_RELAX^time step relaxation factor ( - 1.01) TIME_TOTA^total time for simulation to reach St. st. CEPS^concentration criterion for steady state NR^6 of grids in radial direction NX^6 of grids in axial direction PER_MEM^membrane permeability CONDECS^ECS conductivity VISC_L^viscosity in lumen R_MID^radial grid location half-way in between nodes X_MID^axial grid location half-way in between nodes N_X_OUT^6 of out put data in axial direction N_R_OUT^6 of out put data in radial direction DIFF^protein diffusivity NFIBRE^6 of fibres in reactor PM^molecular weight of protein [g/L] VIRA_I^1st virial coeff. for osmotic relation VIRA_2^2nd virial coeff for osmotic relation RG^gas constant T^temperature [K] TRA_TIME_OUTfield with transient reporting times C_START^average start concentration in ECS IMPLICIT REAL*8 (A-H2O-Z) CALL SET_UP START_TIME = F_TIME0 CALL HEART CPU_TIME = F_TIMEO-START_TIME WRITE(",100) CPU_TIME/60. STOP 100 FORMAT( ' CPU-time: ', F12.4,' [min] ) END SUBROUTINE HEART C This routine calls all other routines. IMPLICIT REAL•8 (A-H2O-Z) PARAMETER( MAXX=20I, MAXR=I I 1 COMMON / B_N / NR, NX, NRM, NXM COMMON / B_DI / RI, R2, R3, XL COMMON / N_OUTP / N_X_OUT, N_ROUT COMMON / B_P2 / P_ECS(MAXX,MAXR) COMMON / CI / T_RELAX, T1ME_TOTA, DT, C_EPS COMMON / CONC / C_ECS(MAXX,MAXR) COMMON / TR_OUT / N_TRANS, TRA_T1ME_OUT(99) N_FILE = 1 TIME_TOTA =0.D0 C_DIFMAX = I.D10 WRITE(20,125) 0, TIME_TOTA CALL UP_DATE_P_OS CALL UP_DATE_P_ECS P_DIFMAX, ITER_P ) WRITE( ,55) ITER_P, P_D1FMAX WRITE(20,55)1TER_P, P_DIFMAX CALL VEL_ECS CALL DATA_OUT ( 0 , TIME_TOTA ) WHILE ( C_DIFMAX .GT. C_EPS ) DO CALL UP_DATE_C ( DT, C_DIFMAX , PE_X_MAX, PE_R_MAX ) C_DIFMAX = C_DIFMAX / DT CALL UP_DATE_P_OS ! up-date the osmotic pressures CALL UP_DATE_P_ECS ( P_DIFMAX,ITER_P ) CALL VEL_ECS^! up_date the velocities in the ECS TIME_TOTA = TIME_TOTA + DT DT = DT • T_RELAX T1ME_H = TIME_TOTA /3600. ! convert to hours IF ( TIME_TOTA .GE. TRA_TIME_OUT( N_F1LE ) ) THEN CC^We reached the time for writing all information WRITE(20,125) N_FILE, TIME_H WRITE(20,130) C_DIFMAX,P_DIFMAXITER_P, PE_X_MAX, PE_R_MAX CALL DATA_OUT ( N_FILE , TIME_H ) IF ( TIME_TOTA .GE. TRA_T1ME_OUT( N_TRANS ) ) THEN CC^we reached the specified max modelling time and we stop! WRITE(' ,200) C_DIFMAX WRITE(20,200) C_DIFMAX RETURN END IF N_FILE N_FILE + 1 END IF WRTTE(.,135)TIME_H, C_D1FMAX, P DIFMAX, ITER_P END WHILE CC We are done the required steady-state cond is obtained CALL DATA_OUT ( N_FILE , TIME_H RETURN • -1--F-1,-.-+-F-F-1-1-4,-1--1--1--1-1-1,+++++++++++++++ +MC++ 55 FORMAT(/3X,'6 of iter. to adjust pressures',16,' tip ',E10.4, ' [Pa]) 100 FORMAT(5X,110,F16.2,1s]) 110 FORMAT(11F7.1) 112 FORMAT(11F73) 125 FORMAT(//' file time:',F9.4,' [h] 130 FORMAT(2X,2(2X,E12.4),2X,17,2(2X,E12.4) ) 135 FORMAT(2X,F10.2,2X,E10.4,2X,E10.4,2X,15,2X,E10.4,2X,E10.4) 200 FORMATCsteady-state was not achieved 1 diffmax C:',E10.4 ) C^WRITE(*, I 12 ) ( (C_ECS(1,1),1=1,NX,N_X_OUT),1=NR, I ,-N R_OUT) C^WRITE(",110 ) ( (P_ECS(1,1),1=1,NX,N_X_OUT),J=NR, I ,-N_R_OUT) END SUBROUTINE UP_DATE_P_ECS ( P_DIFMAX, ITER) • • • The Laplace eqn. that governs the pressure distribution in the ECS space is solved by a line by line method. d2P/dx2 + 1/rd ( r dP/dr) dr = 0 IMPLICIT REAPS (A-H2O-Z) PARAMETER ( MAXX=20I, MAXR=1 I) COMMON / B_N^/ NR, NX, NRM, NXM COMMON / B_PA / P_EPS, NITER P, ALPHAX, ALPHAR COMMON / Bp!^/ PLU(MAXX), P_OS(MAXX) COMMON / B_P2 / P_ECS(MAXX,MAXR) COMMON / PRESS_T / TAX(MAXX), TCX(MAXX), TAR(MAXR), TCR(MAXR), T_F1ELD(MAXX, MAXR) COMMON / C_I^/ T_RELAX, TIME_TOTA, DT, C_EPS COMMON / DUMMY / CONSTANT DIMENSION TB(MAXX), TD(MAXX), P_NEW(MAXX) CC CONSTANT= RI • PER_MEM / ( COND_ECS • A) ALPHAXM = I.DO - ALPHAX ALPHARM = 1 DO - ALPHAR DO ITER = 1, NITER_P P_DIFMAX = 0 DO 1st sweep as r direction DO I = I, NX IF =I + 1 IM =I - I DO = 1, NR TB(.1) = T_FIELD(1,1) TD(1) = P_ECS(IM,11" ( -TAX(I)) 100 + P_ECS(IP„.1) • ( -TCX(I) ) END DO TB(l) = 713(1) - CONSTANT Include BC. at the membrane/ECS interface TD(t) = TD(1) - CONSTANT • ( P_LU(I) + P_OS (D ) CALL TDMA ( TAR, TB, TCR, ID, P_NEW, 1, NR) DO I= 1, NR P_NEW(1) = ALPHARM • P_ECS(1,1), ALPHAR • P_NEW(J) P_DIFMAX = DMAX1(P_DIFMAX, DABS(P_NEW(1)-P_ECS(1,1))) PECS(1,1)= P_NEW(1) END DO END DO C^ 2nd sweep in x direction C---- incorp. the B.C. into the first row =1 1P = .1+ 1 DO! = I, NX TB(1)=T_FIELD(1,1) TD(!) = - CONSTANT • ( P_LU(I) - P_ECS(1,1) + P_OS(I) ) + P_ECS(I,JP) * ( -TCR(I) ) END DO CALL TDMA (TAX, TB, TCX, TD, P_NEW, I, NX ) DO! = I, NX P_NEW(I) = ALPHAXM * P_ECS(1,1) + ALPHAX * P_NEW(I) P_DIFMAX = DMAXI (P_DIFMAX, DABS(P_NEW(1)-P_ECS(1,1))) P_ECS(1,1)= P NEW(I) END DO DO I = 2, NR JP =I + I JM I - I DO! = I, NX TB(I) = T_FIELD(1,1) TD(I) = P_ECS(I,JM) • (-TAR(J) ) + P_ECS(1,1P) • ( -TCR(1) ) END DO CALL TDMA ( TAX, TB, TCX, TD, P_NEW, I, NX ) DO! = I, NX P_NEW(I) = ALPHAXM • P_ECS(1,1) + ALPHAX • P_NEW(I) P_D1FMAX = DMAX1(P_DIFMAX, DABS(P_NEW(1)-P_ECS(1,1))) P_ECS(1,1)= P_NEW(I) END DO END DO CC^now we up-date the lumen pressures CALL LUMEN_PRESS ( P_DIF_L ) CC^check if the requ. ace. is obtained IF ( ( P_DIFMAX + P_DIF_L )^PEPS) RETURN END DO CC We exceeded the maximum iteration and print out an warning WRITE( •,120) T1ME_TOTA/3600., P_DIFMAX, P_DIF_L WRITE(20,120) TIME_TOTA/3600., P_DIFMAX, P_DIF RETURN 110 FORMAT(11F7.1) 120 FORMAT(3)Cpressure diff at t=',F12.3,(h) is stilt,/ . 2X'ECS ',E12.5,/ 2X'Lumen ' ,E12.5 ) END SUBROUTINE LUMEN_PRESS ( P_DIFMAX ) CC ^ This S.R. calculates the new lumen pressures. CC^d2P 1 /dx2 = 16 Lp/R I ^3 • Pt I6Lp/R 1 ^3 ( Pi3 - P2 ) IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=201, MAXR=11 ) COMMON! B_N / NR, NX, NRM, NXM COMMON! B_A / X (MAXX),R (MAXR), D_X(0:MAXX), D_R(0:MAXR) COMMON / B_DI / RI, R2, /13, XL COMMON / B_SP / PER_MEM , COND_ECS, VISC_L COMMON! B_P1 / P_LU(MAXX), P_OS(MA XX) COMMON! B_P2 / P_ECS(MAXX,MAXR) COMMON! PRESS / PIN, P_OUT DIMENSION P_NEW(MAXX), A(MAXX), B(MAXX), C(MAXX), D(MAXX) C^ now we calculate the new axial pressures in the lumen with a C^ Finite Difference Method. P_DIFMAX = O.DO DUMMY= 16.D0 • PER MEM / R I **3.D0 A(1)= O.DO B(I)= LDO C(I) = 0.130 0(l) = P_1N ! pressure at the fibre lumen inlet DO I = 2, NXM R_DX1M = I .00 / D_X(I-1) R_DXI = 7130 / D_X(I) FISUMDX = ( D_X(I-1) + D_X(I) ) • 0.500 A(I) =^R_DXIM B(I) = - ( R_DXI + R_DX1M + DUMMY " FISUMDX ) C(I) = R_DX I D(I) = DUMMY • HSUMDX • P_OSM-P_ECS(1,1) END DO A(NX)= O.DO B(NX) 1.00 C(NX)= 0.00 D(NX)= P_OUT ! pressure at the fibre lumen outlet CC^solve the tridiagonal matrix and up-date the lumen pressures CALL TDMA ( A, B, C, D, PNEW, I, NX ) DO I = 2, NXM P_DIFMAX = DMAXI( P_DIFMAX, DABS( P_NEW(I)-P_LU(I))) P_LU(I) = P_NEW(I) ! new lumen pressure END DO RETURN END SUBROUTINE DATA OUT ( NUM, TRANS:TIME ) C ^ This Function writes the data to a file IMPLICIT REAL*8 (A-H2O-Z) CHARACTER*3 INT_CHAR, NUMBER PARAMETER ( MAXX=20I, MAXR=Il ) COMMON / B_VE_ME / VXR_LU(MAXX), UXR_LU(MAXX) COMMON / B_N /NR, NX, NRM, NXM COMMON! B_P I / P_LU(MAXX), P_OS(MAXX) COMMON! B_P2 / P_ECS(MAXX,MAXR) COMMON! B_A / X (MAXX),R JMAXR), D_X(0:MAXX), D_R(0:MAXR) COMMON / CONC / C_ECS(MAX)C,MAXR) COMMON! B VE ECS/ VEL(MAXX,MAXR,4) COMMON! FIALF-DR / X MID(MAXX+ I), R_MID(MAXR+ I) COMMON / B_SP / PER MEM , COND_ECS, VISC L COMMON / BDI / RI, R2, R3, XL COMMON! AVER / C_AVER(MAXX) COMMON / N_OUTP / N_X_OUT, N_R_OUT NUMBER = INT_CHAR( NUM) 101 = 11 102 = 12 CC^calculate a scaling factor for the velocities inin/h FACTOR = 1000 • 3600 OPEN( UNIT= 101, FILE = 'D2_DAT_ANNUMBER ) OPEN( UNIT=IO2, FILE = D2_DAT_BY/NUMBER ) CC^Calculate a mass balance in the ECS CALL C_AVERAGE CC^Calculate the velocities in the lumen CALL LUMEN_VEL CC^First we write the averaged values to separate file DO I = 1, NX, N_X_OUT WRITE(I01,100) Mr 100., C_AVER(I), UXR_LU(I), VXR_LU(1)*FACTOR, F_PI( C_AVER(1) ), P_LU(I), F_R_AVRAGE( P_ECS,I ) END DO 100 FORMAT ( 2(F8.4,2X), 2(E14.6,2X), 3(F I2.4,2X) ) DO 3=1, NR, N_R_OUT JM =1 - I DO I = 1, NX, x_x_our IM = I - I WRITE(IO2,101) X(1)•1132, R(J)•1136, C_ECS(I,1), P_ECS(1,1), CC^interpolate the axial velocity in the ECS at the grid-points (D_X(1)*VEL(1,1,2)+D_X(IM)*VEL(1,1,1))/(D_X(1)+D_X(IM)) *FACTOR, CC^interpolate the radial velocity in the ECS at the grid-points (D_R(.0'VEL(1,1,4)+D_R(IM)•VEL(t,1,3))/(D_R(1)+D_ROMD *FACTOR END DO END DO WRITE(I01,110) TRANS_TIME ! time-stamp on each file WRRE(IO2,110) TRANS_T1ME ! time-stamp on each file CLOSE( UNIT=101 ) CLOSE( UNIT=IO2 ) RETURN 101 FORMAT ( 3(F8.4,2X), F12.4,2X, 2(E14.6,2X)) 110 FORMAT ( /F14.3) END SUBROUTINE UP_DATE_C ( TIME_STEP, C_DIFMAX, PE_X_MAX, PE_R_MAX ) CC This routines solves the 2D diffusion convection eq in the EC S CC The PDE eq is discretized by a finite difference approach CC and to ensure realistic solutions a "power-law" scheme with 101 CC 1st order up-winding is included. CC At the first half-time step the radial concentration is calc and CC in the next half-time step the axial concentration is updated. This CC method is also referred to as the ADI-Method, for further details see Anderson CC eta!. 1980. IMPLICIT REAL*8 (A-H2O-Z) PARAMETER( MAXX=20I, MAXR=I I ) COMMON! B_N / NR, NX, NRM, NXM COMMON / BA / X (MAXX),R (MAXR), D_X(0:MAXX), D_R(0:MAXR) COMMON! B_DI / RI, R2, R3, XL COMMON / HALF_GR / X_M1D(MAXX+ I), R_MID(MAXR+ I) COMMON! CON / A_R(MAXR), A_X(MAXX) COMMON! B VE ECS / VEL(MAXX,MAXR,4) COMMON! CONC- / C ECS(MAXX,MAXR) COMMON / DIFFUS / DiF F,NFIBRE DIMENSION A(MAXX), B(MAXX), C(MAXX), D(MAXX), CNEW(MAXX), . PE_R(MAXX,MAXR), PE_X(MAXX,MAXR), CSTAR(MAXX,MAXR), AV(MA30), AN2(MAXX,MAXR), AS2(MAXX,MAXR),AW2(MAXX,MAXR),AE2(MAXX,MAXR) C_D1FMAX = 0.D0 PE_R_MAX = 0.D0 PE_X_MAX = 0.D0 DT = 2.0D0 / TIME_STEP^! set rezipr. halftime-step for AD! C set-up the Peclet numbers in axial PE_X and radial PE_R direction DO 1 -= I, NR DO 1= I, NXM PE = VEL(1,1,1)* D_X(I)/ DIFF PE_X(1,1)= F_POWER( PE ) • DIFF / D_X(I) calculate the highest local axial Pe number PE X MAX = DMAXI(PE_X_MAX, FE) END- D-0 END DO DO I = 1, NX DO J = 1, NRM PE = VEL(1,1,3)* D_R(1) / DIFF PE_R(1,1)= F_POWER( FE) * DIFF/ D_R(1) calculate the highest local radial Pe number PE_R_MAX = DMAX1( PE_R_MAX, FE) END DO END DO CC Now we calculate some repeated constants DO I = I, NX 1M =I - I DO1 = 1, NR JM =1 - I JP I + I AN2(1,J) = ( PE_R(1,1 ) + DMAX I (-VEL(1,1,3),O.D0) ) • R_MID(JP)*A_R(J) A52(1,1) = ( PE_R(I,JM) + DMAX1( VEL(1,1,4),0.D0) ) • R_MID(1) • A_R(.1) AW2(1,1)= ( PE_X(IM,J)+ DMAXI( VEL(1,1,2),0.D0) ) •A_X(I) AE2(I,J) = ( PE_X(I,1 ) + DMAX I (-VEL(I,1,1),O.D0) ) • A_X(1) END DO CC^This array spec. the constants at the membrane/ECS interface AV(I) VEL(1,1,4) * A_R(I) • R2 END DO CC Now we start with the ADI-Method. • ^ CC At the 1st half-time step we sweep over all x grids to cats ^ CC^the new concentrations in radial direction 00 1= = I, NX IM = I - I IF = I + 1 =I A(J) = 0.00 C(J) = -AN2(I,J) 13(1) = DT - C(J) - A(J) + AVM D(1)=^AW2(I,J) * C_ECS(IM,J) + ( DT-AE2(1,1)-AW2(1,1) ) * C_ECS(1,1) AE2(1,1) * C_ECS(IP,1) DO 1 = 2, NRM JM =1-1 JP = 1+1 A(1)= - AS2(I,1) CM= - AN2(1,J) B(J) = DT - C(i) - A(1) D(J)= AW2(1,1)^" C_ECSHM,11 + ( DT- AE2(1,1) - AW2(1,1) C_ECS(1,1) +^AE2(1,1)^• C_ECS(IP,J) END DO =NR 1M = NR - I A(NR) = - AS2(1,1) C(NR) = 0.D0 B(1) = Dl- C(1) - AU) D(NR) =^AW2(1,1) • C_ECS(IM.11 +( DT - AE2(I,1) - AW2(l,F) ) • C_ECS(1,1) +^AE2(I,1)^• C_ECS(IP,1) CALL TDMA ( A, B, C. D, CNEW, I , NR) DO 1= I, NR CSTAR(1,1) = CNEW(1) END DO END DO CC Note: the intermediate concentration values CSTART are meaning less CC^and only used for the next time step CC Now in the 2nd half-time step we sweep over all r grids to cats. ^ CC the new concentrations in axial direction. =1 ! first row JP = + I DO I = I, NX A(1)= -AW2(1,1) -AE2(I,J) • B(1)= DT - C(I) - A(1) • (DT - AN2(I,J) - AVM ) * CSTAR(1,1) AN20,1)^• CSTAR(1,1P) + END DO CALL TDMA ( A, B, C, D, CNEW, I, NX) DO I = I, NX C_DIFMAX = DMAXI(C_DIFMAX,DABS(CNEW(I)-C_ECS(1,1)) ) C_ECS(I,H= CNEW(I) END DO DO J = 2, NR JM =1 - I JP =1+ I DO I = I, NX A(I) = -AW2(1,1) C(I)= -AE2(I,J) BM= DT - C(I) - A(I) D(I) =^AS2(I,J) CSTAR(1,1M) + ( DT - AN2(1,1) - AS2(1,1) ) * CSTAR(1,11 • AN2(I,J)^CSTAR(1,1P) END DO CALL TDMA ( A, B, C, D, CNEW, I, NX) DO I = I, NX C_DIFMAX = DMAX1(C_DIFMAX,DABS(CNEW(1)C_ECS(1,1)) ) C_ECS(1,1) = CNEW(I) END DO END DO RETURN END REAL*8 FUNCTION F_POWER( Pc) C This function applies the power law scheme for the finite C difference formula. Other schemes can be used by substituting C appropriate functions (see Pantankar (1980) for details]. IMPLICIT REAL•8(A-H2O-Z) F_POWER = DMAX1( 0.D0 , I.D0 - 0.1 DO *I DABS(Pe)^) RETURN END SUBROUTINE SET_UP IMPLICIT REAL•8 (A-14,0-Z) PARAMETER( MAXX=20I, MAXR=Il ) COMMON! B_P1 / PI COMMON / B_A / X_(MAXX'),R (MAXR), D_X(0:MAXX), D_R(0:MAXR) COMMON! B_PA P_EPS, NITER P. ALPHAX, ALPHAR COMMON / B_PI / P_LU(MAXX), P_OS(NIAXX) COMMON! B_P2 / P_ECS(MAXX,MAXR) COMMON! C_I / T_RELAX, TIME_TOTA, DT, C_EPS COMMON! B_N / NR, NX, NRM, NXM COMMON / B_DI / RI, R2, R3, XL COMMON! DUMMY / CONSTANT COMMON / B_SP / PER_MEM, COND_ECS, VISC_L COMMON / PRESS _T / TAX(MAXX), TCX(MAXX), TAR(MAXR), TCR(MAXR), T_FIELD(MAXX, MAXR) COMMON / HALF OR / X_MID(MAXX+ 1 ), R_MID(MAXR+ 1 ) COMMON / N_OUTP / N_X_OUT, N_R_OUT COMMON! CONC / C_ECS(MAXX,MAXR) COMMON / DIFFUS / DIFF, NFIBRE COMMON / OSMOTI /PM, T, VIRA_I, VIRA_2, RG, OS CON COMMON / CON / A_R(MAXR), A_ X(MAXX) COMMON ! TR_OUT % N_TRANS, TRATI ME _our(99) COMMON / C_BEGIN / C_START 102 COMMON / PRESS / PIN, POUF OPEN( UNIT-10, FILE = 'HF2D_IN.PUP) OPEN( UNIT-20, FILE ='D2_DARCY.INF) READ(10,*) NX !8 of grids in X-direction READ(I 0,*) NR !8 of grids in R-direction READ(I0,*) N_X_OUT ! 14 of grids in X-direction READ(I 0,*) FLOW ! flow rate in the reactor [mL/min] READ(I0,*) P_IN ! Inlet pressure [Pa] READ(10,*) P_OUT ! Outlet pressure [Pa] READ(I0,*) XL ! [m] length of the capillary READ(10,*) RI^! [m] inner lumen radius READ(I0,*) R2^[rn1 outer membrane radius READ(I0,*) R3 ! [m] Krogh radius FtEAD(I0,*) RSHELL ! [m] inner shell radius of the HFBR READ(I 0,*) NFIBRE ! [-] 8 of fibres in the reactor READ(I0,*) PER_MEM ! [m] membrane permeability READ(I0,*) COND_ECS! [mY] ECS conductivity READ(I0,*) C_START ! BSA start concentration in the ECS READ(10,*) DIFF ! BSA diffusivity READ(I0,*) PM ! molecular weight of the protein READ(10,•)TC ! temperature [Co] READ(I0,*) G_FX ! Grid-factor for sym. grid READ(I0,*)0_FR ! Grid-factor for unsynt. grid READ(I 0,*) P_EPS ! Accuracy for P_ECS in ADI READ(I0,*) ALPHAX ! SOR factor for pressure in x-direction READ(I0,*) ALPHAR ! SOR factor for pressure in R-direction READ(I0,*) NITER_P ! max. 8 of iterations for pressures READ(I0,*)C_EPS ! steady-state acc. for the concentrations READ(I0,*) DT ! Time-step for the cliff cony. eq. READ(10,*)T_RELAX ! time-relaxation factor for the dill'. cony. eq. T^= 273.15 + TC ! convert to Kelvin V1SC_L = F_VISC_H20( TC ) CC^viral coeff for osm. pressure relationship VIRA _I = 0.010473D0 VIRA_2 = I.73743D-5 RG^= 8.3141500^! GAS CONSTANT.Umoles-K OS_CON = RG • T / (PM/1000.) DO I = I, 99 READ(l 0,•,ERR=18) TRA_TIME_OUT(I) END DO 18^N_TRANS = I - I CLOSE( UNIT=I0 ) C^convert to seconds ! DO I = 1, N_TRANS TRA_TIME_OUT(I) = TRA_TIME_OUT(I) • 3600.00 END DO IF ( NX .GT. MAXX) NX=MAXX IF ( NR .GT. MAXR) NR=N1AXR CC^Set the relaxation parameter for the SOR in line by line if they are CC^not appropriate ! PI = DATAN(I.D0) • 4.00 IF ( ALPHAX .LE. 0.00 .0R. ALPHAX .GT. 200) THEN ALPHAX=2./(1.+DSQRT(I -(1. DO - PI**2./(2*NX**2))**2)) END IF IF ( ALPHAR .LE. 0.00 OR. ALPHAR .GT. 2.00) THEN ALPHAR = 2./(1.+DSQRT(1-(1.D0-P1**2./(2*NR••2))**2) ) END IF WRITE(*,32) ALPHAX, ALPHAR 32 FORMAT(I0X,'Relaxation Parameters:',/, 5x, 'alpha-x: ',F7.4,' alpha-r : ',F7.4) CC^Check if the Krogh Radius is given !! IF ( R3 .LE. 0.00) THEN R3 = RSHELL / NFIBRE**0.5 WRITE(*,*) 'Krogh radius:', R3 END IF IF ( FLOW .GT. 0.000) THEN P_OUT = 0.00 QF = FLOW*I.D-6/60/NFIBRE ! calculate flow in a single fibre P_1N = QP`XL*8*VISC_L/(P1*R I **4.) END IF WRITE(,100) P_IN, P_OUT 100^FORMAT(5X,'p-in:',F10.2,5X;p-out:',F I 0.2,' [Pap FLOW = P_IN/(1.D-6/60/NFIBRE•XL*8•VISC_D(Pl•RI".4.)) WRITE(*,110) FLOW 110^FORMAT(5X,'Q-HFBR:',F I 0.2; [ml/min]') DO I = I, NX PLUM = F_YLIN( 0.00, XL, P_IN, P_OUT, X_(I) ) P OS(I) = F_P1( C_START ) D-0 1= 1, NR P ECS(I,1)= DP ECS(1,1) = C_START END DO END DO PLUM = P IN P_LU(N)G =POUT C---- set-up coeff for the pressure eq. DO I = 1, NX 1M^= I-1 X_IN = ( D_X(I) + D X(IM) ) • 0.5D0 ! integration from w to e X_W = X_(I) - D_X-(IM) 0.5D0 X_MID(I) = X W A X(I) = 1.130—/X IN IF—( 1. EQ. I ) THEN TAX(1)= 0.00 ELSE TAX(I)= 1.00 / ( D_X(1M) • X_IN END IF IF ( I. EQ. NX ) THEN TCX(I) = 0.00 ELSE TC X(I) = I .D0 / ( D_X(I) • X_IN ) END IF END DO X MID(NX+ I) = X_(NX) DO 3 = I, NR JM = 1-1 R_S = R_(.1) - D ROM) • 0.500 R_N^R_(J) + 13_R(.1 ) • 0.5D0 R M1D(J) = R S DUMMY =W2. - R_S2.)* 0.500 A_R(1) = 1./DUMMY A^= DUMMY^! integration area IF ( J. EQ. I ) THEN TAR(J) = 0.00 CONSTANT — RI • PER_MEM / ( COND_ECS • A) ELSE TAR(J) = 1.00 / (DLOG( R(1)/(R(1)-D_R(1M)) ) • A) END IF IF ( 1. EQ. NR ) THEN TCR(J) = 0.00 ELSE TCR(.0 = 1.00 / (DLOG( (R_(.1)+ D_R(1))/R (1)) " A) END IF END DO R MID(NR+I) = R_(NR) DO I = 1, NX DO I = 1, NR T FIELD(1,1) = - TAR(J) -TCR(1) - TAX(I) - TCX(I) END DO END DO CALL DATA SHOW( RSHELL, FLOW, GF_X, GF_R ) RETURN END REAL*8 FUNCTION F_VISC_H20( TC ) C This function calculates the viscosity of H2O at the given Temp TC in Celsius. C input TC : oC C output vise : Pus IMPLICIT REAL*8(A-H2O-Z) IF (TC .GT. 99.00 .0R. TC .LE. ODO ) THEN WRITE() 'EXCEEDING RANGE IN VISCOS EQ. II!' END IF F_VISC_H20 = 10.1:10"*( (1.3272*(20-TC)-0.001053*(TC-20)••2 ) / (TC+I05.) ) • 1.0021/100/10 RETURN END REAL*8 FUNCTION F_VISC_BSA( C, VISC TIME_TOTA = 0.00 P_TIME_TOTA = 0.00 NXM = NX - I NRM = NR - 1 N_X_OUT = NXM / (N_X_OUT-1) N_R_OUT = 1 CC ^ Setup the grids in both directions CALL USYM_GRID ( 0.DO, XL, NX, G_FX, X_, D_X ) CALL USYM_GRID ( R2 , R3, NR, G_FR, R_, DR) CC ^ Initialize a linear pressure drop along the lumen-side! DP = (P_IN - P_OUT)*0.5D0 I average pressure in the ECS C This function relates the viscosity to the BSA concentration. C The viscosity of water must be also specified C input C BSA g,'L VISC H20 : Pa s C output visc : Pus IMPLICIT REAL•8(A-H2O-Z) V_INT = 3.65D-3 V_K = 2.08500 C2 = VINT * C F_VISC_BSA = VISC*( I a C2 * ( I V_K *C2 ) ) RETURN END 103 100 FORMAT(5X;mass-balance ',F8 4,' lost.',F10.6; Ig/LI'l RETURN END SUBROUTINE VEL_ECS CC This routine up-dates the velocities in the ECS. The velocities CC^are located half-way between the grids CC IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=201, MAXR=I I) COMMON / B_N / NR, NX, NRM, NXM COMMON / B_A / X JMAXX), R_(MAXR), D_X(O.MAXX), D_R(0:MAXR) COMMON / B_DI / RI, R2, R3, XL COMMON! B_SP / PER_MEM, COND_ECS, VISC_L COMMON! B_PI / P_LU(MAXX), P_OS(MAXX) COMMON! B_P2 / P_ECS(MAXX,MAXR) COMMON! B_VE_ECS/ VEL(MAX)çMAXR,4) COMMON / HALF_OR / X_MID(MAXX+ I), R_MID(MAXR+1) COMMON! CONC / C_ECS(MAXX,MAXR) DUMMY = -COND_ECS/VISC_L CC^calculate velocities in r-direction V(x,r) n and s DO .1= 2, NRM, 2 JP = 1+1 JM =1-1 DO I = I, NX VEL(1,1,4) = DUMMY * (P_ECS(1,J)-P_ECS(1,1M))/D_R(.1N1) VEL(1,1,3) = DUMMY • (P_ECS(1,JP)-P_ECS(1,1))/D_R(1) VEL(1,1M,3)= VEL(1,1,4) VEL(1,111,4)= VEL(1,1,3) END DO END DO CC calculate velocities in x-direction U(x,r) e and w DO I= 2, NXM, 2 IP = 1+1 1M = I-I DOS = I, NR VEL(I,J,2) = DUMMY * (P_ECS(1,J)-P_ECS(1M,J))/D_X(IM) VEL(1,1, I ) = DUMMY • (P_ECS(lP,J)-P_ECS(1,1))/D_X(1) VEL(tM,1,1)= VEL(1,1,2) VEL(IP,1,2)= VEL(1,1,1) END DO END DO DUMMY = RI • PER MEM / ( VISC_L * R2) DO! = I, NX CC^velocity on the grid-points at the membrane out-side R2 VEL(I, 1,4) = DUMMY • ( P_LU(1)-PECS(1,1)+P_OS(1) ) CC^VEL(1,NR,3)= 0.D0 ! n at R3 END DO DO .1= I, NR VEL(I, 1,2) = 0.0D0 ! at X = 0 W VEL(NX.,1,1)= 0.0D0 ! at X = L E END DO RETURN END ******* ****************^............... **•**■• frt*• • SUBROUTINE C_AVERAGE C C This S. R. calculates the average concentration in the radial directions. IMPLICIT REAL"8 (A-H2O-Z) PARAMETER ( MAXX=201, MAXR=I t) COMMON! B_A / X JMAXX),R JMAXR), D_X(0:MAXX), D_R(0:MAXR) COMMON! B_N / NR, NX, NRM, NXM COMMON! CONC / C_ECS(MAXX,MAXR) COMMON! AVER! C_AVER(MAXX) COMMON / C_BEG1N / C_START CC Check for negative concentrations DO I = I, NX DO .1= I, NR IF ( C_ECS(1,5). LT. 0.1)0) C_ECS(1,1) = 0.1)0 END DO END DO CC Calculate the average radial value on each axial grid location DO I= I, NX C_AVER(I) = F_R_AVRAGE( C_ECS, 11 END DO CC Now we calculate a mass-balance over the reactor CC^First we calculate a spline over the radial averaged concentrations CALL SPLINE ( X_ , C_AVER, NX, 2) CC CC now we integrate the area under this spline and divide by the reactor length to obtain the total average concentration in the ECS CAV = FS_INTEG ( X_( I ), X (NX) ) / X JN WR1TE(20,100) CAV , C_START-CA V WRITE( ',100) CA V , C_START-CAV REAL'8 FUNCTION FR _AVRAGE FIELD, 1 CC This function returns a radial average value for any spec x location. IMPLICIT REAL•8 (A-H2O-Z) PARAMETER MAXX=201, MAXR=11 COMMON / B_A / X JMAXX),R JMAXR), D_X(0:MAXX), D_R(0 MAXR) COMMON! B_N / NR, NX, NRM, NXM DIMENSION FIELD ( MAXX, MAXR ) SUM =0.1)0 C23 = 2.D0/3.D0 DO I = 2, NR JM =1 - 1 SUM = (^D_R(JM) FIELD(1,1M) + 2.1)0 R_(JM) • FIELD(1,1M) + C23 D_R(JM) ( FIELD(1,1)-FIELD(1,JM) ) + R J1M) • ( FIELD(1,J)-FIELD(1,1M) ) ) • D_R(JM) + SUM END DO CC now we divide this by the annular cross-section to obtain the ^ CC required average value F_R_AVRAGE = SUM! ( R (NR)"•2 - R_(I)"2 ) RETURN END SUBROUTINE DATA_SHOW ( RSHELL,FLOW,GF_X,GF_R ) CC This routine writes the important modelling parameters to an file IMPLICIT REAL*8 (A-H2O-Z) PARAMETER( MAXX=201, MAXR=I I ) COMMON! B_PI / PI COMMON / B_N I NR, NX, NRM, NXM COMMON / N_OUTP / N_X_OUT, N_R_OUT COMMON! B_DI / RI, R2, R3, XL COMMON! B_SP / PER MEM, COND_ECS, VISC_L COMMON! DIFFUS! DIFF, NFIBRE COMMON! OSMOT I / PM, T, V1RA_ I, VIRA_2, RG, OS CON COMMON / TROUT / N_TRANS, TRA_TIME_OUT(99) COMMON! C_BEGIN/ C_START COMMON! PRESS / PIN, P_OUT COMMON / B_PA / P_EPS, NITER_P, ALPHAX, ALPHAR COMMON! BA / X (MAXX),R JMAXR), D_X(0:MAXX'), D_R(0:MAXR) COMMON! C_1 / T_RELAX, TIME_TOTA, DT, CEPS 10 = 33 OPEN( UNIT=I0, FILE = 'D2_DARCY.DAT) WRITE(10,001) NX,NR 001 FORMAT(6X,14,' x-grids ',14,' r-grids) WRITE(I0,003) PIN, P_OUT 003^FORMAT(2X;Pressure: inlet ^:',F10.2; War/ . 2X,'^outlet^1',F10.2; 'Par) WRITE(I0,005) FLOW 005 FORMAT(2X;Vol. flow rate in HFBR^:',F10 4,' [ml/minr) WRITE(10,007) XL*100. 007 FORMAT(2X,'Fibre length^:',F10.4,. (cm)) WRITE(10,009) RI " I D6, R2•1D6 009 FORMAT(2X;Fibre radii- inner^:',F10.2; ()am)'! outer^:',F10.2; ltrnir) WRITE(I0,01 I) R3.1D6 011 FORMAT(2X;Krogh radius ^:',F10.2,' 4un(') WRITE(10,013) RSHELL^100. 013 FORMAT(2X;Inner cardrige radius ^.',F I 0.4; (emi(') WRITE(10,015) NFIBRE 015 FORMAT(2X;# of fibres in reactor^:',110) WRITE(I0,016) NFIBRE*R2•XL•PI•2 016 FORMAT(2X,'outer HFBR membrane area^ ',F10.4,1m2r) WRITE(I0,017) VISC_L 017 FORMAT(2X,'Viscosity in lumen ^ ',E10.4; [Pa-s(') WRITE(I0,019) PER MEM 019 FORMAT(2X;Meurb rane permeability^.',E10 4,' (inn WRITE(I0,020) PER MEM^RI • DLOG )R2/12.1) 020 FORMAT(2X;Menrbrane conductivity^:',E10.4; (11121') WRITE(10,021 ) COND_ECS 021 FORMAT(2 X;ECS-Conductivity^/.E10 A,' [m2(') WRITE(I0,023) 023 FORMAT(/,10X,' Protein specifications in ECS: WRTTE(10,025) C_START 025 FORMAT(2X,'Start Concentration^.',F10.4,' [g/L(') WRITE(10,027) DIFF 027 FORMAT(2X;Diffusivny^:',E I 0.4,' (m2/sr) WRITE(10,029) PM 104 029 FORMAT(2X,'Molecular weight ^:910.2, [Dap WRITE(10,031) F_PI(C_START) 031 FORMAT(2X,'Osm press. at start^:',910.4,' (Pa]) WRITE(10,032) F_V1SC_BSA(C_START, V1SC_L) 032 FORMAT(2X,'Viscosity of protein so! . :',E10.4,' [Pa-s]') WRITE(I0,033) T -273.15 033 FORMAT(2X,Temperature^:',F10.2,' [oC]',/) WRITE(I0,037) PEPS 037 FORMAT(2X,'Pressure tolerance^WA,' [Pa]') WRITE(10,039) NITER_P 039 FORMAT(2X,'Max iter. to reach the toler_ :',I10) WRITE(I0,041) ALPHAX, ALPHAR 041 FORMAT(2X,'axial relaxation parameter :',F10.4,/ .^2X;radial^:',F10.4) WRITE(10,043) CEPS 043 FORMAT(2X,'Steady-state for conc. ^:',E I 0.4,' [g/L]') WRITE(10,045) DT, T_RELAX 045 FORMAT(2X„'inital time step^:',F10.4,/ aced. factor:^:910.8) CC This Sub. R. up-dates the osmotic pressures along the membrane/ECS CC interface. IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=20I, MAXR=11 ) COMMON / B_N / NR, NX, NRM, NXM COMMON / B_P / P_LU(MAXX), P_OS(MAXX) COMMON / CONC / C_ECS(MAXX,MAXR) DON I,NX P_OS(I) = F_P1( C_ECS(1,1) ) END DO END REAL 8 FUNCTION F_PI( C_IN ) IMPLICIT REAL*8 (A-H2O-Z) COMMON / OSMOT1 / PM, T, VIFtA_ I, VIRA_2, RG, OS_CON F_PI = OS_CON • C_IN * (I DO + C_IN (VIRA_I+C_IN*VIRA_2) ) RETURN END WRITE(10,049) 049 FORMAT(2X,//,'**• Transient times for print-out *** ) DO! = 1, N_TRANS TR =TRA_TIME_OUT(1) WRITE(10,055)111, TR/60. END DO 055 FORMAT(5X, F10.2, ' [5] ',F10.4,' [min] ) WRITE(I0,057) GF_X 057 FORMAT(2X//,'*** grid loc./spacing in x-dir. [cm] **• ',F10.6) DOI= I, NX WRITE(I0,059) I, X(I)*100., D_X(1).100. END DO 059 FORMAT(2X,I4,2(5X, F10.6) ) WRITE(10,061) GF_R 061 FORMAT(2X1/,'••• grid locations in r-direction [mml ^' .^,F10.6) DO I = I, NR WRITE(10,059) J,^D_R(1)•1.D6 END DO 063 FORMAT(2X,142(5X, F10.6) ) CLOSE( UNIT=10 ) RETURN END ..... ...... *..***** ****** *...5•1•**.**** ********* REAL"8 FUNCTION F_YLIN( XI, X2, Y1, Y2, XLIN ) CC This Function calc. a linear interpolation CC between two points. IMPLICIT REAL*8 (A-H2O-Z) F_YLIN = ( Y2 - Y1 X2 - XI) • ( XL1N - XI) + YI RETURN END CHARACTER•3 FUNCTION INT_CHAR ( I ) **** SUBROUTINE LUMEN_VEL CC calculate lumen velocities IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=201, MAXR=I I) COMMON / BPI / PI COMMON / B_N / NR, NX, NRM, NXM COMMON! B_A / X_(MAXX) ,R_(MAXR), D_X(0:MAXX), D_R(0:MAXR) COMMON! B_D1 / RI, R2,11.3, XL COMMON / B_SP / PER_MEM , COND_ECS, V1SC_L COMMON! B_PI / P_LU(MAXX), P OS(MAXX) COMMON / B_P2 / P_ECS(MAXX,MAXR) COMMON / DIFFUS / DIFF, NFIBRE COMMON / B_VE_ME/ VXR_LU(MAXX), UXR_LU(MAXX) DIMENSION VDUMMY(MAXX) C ^ UXR_LU: contains the average axial velocities in the lumen C^ VXR_LU: contains the radial velocities at RI (inner membrane) C^ calculate the average lumen axial velocities DUMMY= - R1**2.D0 /(8.*VISC_L) DO I= 2, NXM IF = I + I IM =I - 1 SLOPE = (P_LU(IP)-P_LU(IM))/(D_X(1)+D_X(IM)) UXR_LU(I) = SLOPE • DUMMY END DO UXR_LU(1) = DUMMY * (P_LU(2)-P_LU(I) )/ D_X(I) UXR_LU(N X) = DUMMY (P_LU(NX)-P_LU(NX-1)) / D_X(N X-1 ) DUMMY = PER_MEM / VISC_L DO I = 1, NX VXR_LU(1) = DUMMY * ( PLUM - P_ECS(1,I) + P_OS(I) ) VDUMMY(I) = DABS( VXR_LU(1) ) END DO Q_U_IN = UXR_LU(1) * PI • RI • 2. ! Bow rate in the lumen C^calculate the leakage flow into the lumen CALL SPLINE ( X_, VDUMMY, NX, 2) O_V_LEK = 2.DO*Pl•RI * FS_INTEG ( X(1), X_(79X) )/ 2.00 FAC = 1.136`60.*NFIBRE WRITE(20,55) Q_U_IN•FAC, O_V_LEK*FAC WRITE('' ,55) Q_U_IN•FAC, O_V_LEK•FAC RETURN 55^FORMATOX,'Q Lu : ',F10.4,' Q Leak : ',F10.4,' [mUminr) END SUBROUTINE UP_DATE_P_OS This Function changes an Integer to an char, string (0-99) IMPLICIT REAL*8 (A-H2O-Z) C HARACTE R*1 A I , A2, A3 ITEL =1110 Al^=0' A2 = CHAR( ITEL + 48 ) A3 = CHAR( I - ITEL•I0 + 48 ) 1NT_CHAR = A I //A2//A3 RETURN END SUBROUTINE USYM_GRID ( XO, XF, N, G_F, X, DX) C^This S.R. establishes a non-uniform grid_ C^For positive G_F values more grids are on the left side. C For negative G_F values more grids are on the right side_ C^0 < GF_A < PI/2 IMPLICIT REAL*8 (A-H2O-Z) DIMENSION X(N), DX(0:N) PIN = 2.D0 • DATAN(I .D0) IF ( DABS(G_F) .GT. PIN) G_F=0.D0 NM = N - 1 IF (0_F .EQ. 0.00) THEN DXC =(XF-X0)/NM X(1)= XO DO I = 2, N X(1) = X(1-1) + DXC DX(I-1)= DXC END DO ELSE G = DABS(G_F) DO I = I, N IM = I - 1 X(I) = XO + (XF- X0) DTAN( 0• IM/NM )/DTAN(G) IF (I NE. I) DX(IM)=X(I)-X(IM) END DO IF ( G_F^0.D0 ) THEN DOI = I, N DX(I) = X(N-1+1) - X(N-I) END DO X(I) = XO DO 1= 2, N X(1) = X(I-I) + DX(I-1) END DO END IF END IF DX(0) = 0.D0 DX(N) = 0.D0 RETURN END SUBROUTINE SYM _GRID ( XO, XF, N, G_F, X, DX) 105 IMPLICIT REAL*8 (A-H2O-Z) C This S.R. establishes a symmetric non-uniform grid. C^0 < GF_A < PI/2 C^Note GF_A = 0.D0 leads to a uniform grid DIMENSION X(N), DX(0:N) =2.D0 • DATAN(I .D0) IF ( DABS(G_F) .GT. NM = N - I IF ( G_F .EQ. O.DO ) THEN DXC = (XF-X0)/NM X0 DO I = 2, N X(I) = X(I-1) + DXC DX(1-1)= DXC END DO ELSE =DABS(G_F) DO I = 1, N IM =I - I X(I) = X0 + (XF-X0) • 0.5D0 • (I .DO+DSIN( G• (2.D0 * IWNM-I.D0 )) /DS1N(G) ) IF ( I NE. I ) DX(IMX(1)-X(IM) END DO END IF DX(0) = O.DO DX(N) = 0.D0 RETURN END SUBROUTINE TDMA ( A, B, C, D, X, NI, N2) CC This S.R. solves a nidiagonal system of equations. In this CC routine the start and end of the eq. in the matrices needs to CC be defined. IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=20I) DIMENSION A(N2), 8(N2), C(N2), D(N2), X(N2), P(MAXX), Q(MAXX) P(NI)= -C(b11)/ B(N1) Q(N1) D(NI)/B(NI) DO I =N1+1, N2 IM "I-I DEN = A(I) • P(IM)+ B(I) P(I) = -C(I)/ DEN Q(1) = (D(1)-A(1)*Q(1M))/DEN END DO X(N2)= Q(N2) DO 1 = N2-1,N1, -1 X(I) = P(1) " X(1+1) + Q(I) END DO RETURN END REAL*8 FUNCTION F_TIME0 *********** ****************** ********* **..*••■■.**•..***••■ C This function returns the time in seconds after midnight INTEGER*2 hour, minute, second, hundredth, YEAR, MONTH, DAY CALL GE I (1M( hour, minute, second, hundredth ) CALL GETDAT (YEAR, MONTH, DAY) WRITE(20,100) DAY,MONTITYEAR, HOUR,MINUTE,SECOND 100 FORMAT (25C,/,'date: ^ ',3(13'1) F_TIME = ((DBLE( hour ) • 3600.0) + (DBLE( minute) • 60.0) + DBLE( second) + (DBLE( hundredth )/ 100.0)) END SUBROUTINE SPLINE ( X, Y, N, IBC) This subroutine calculates the Q,R,S values for the spline interpolation. The boundary for each end can be selected ot of the three types 1:natural splines 2:clamped splines 3: fitted splines IMPLICIT REAL•8(A-H2O-Z) PARAMETER( MAXX=20I ) COMMON / BLKA / XX(MAXX),YY(MAXX), NN COMMON / BLKB / Q(MAXX), R(MAXX),S(MAXX) DIMENSION X(N), Y(N), H(MAXX) DIMENSION A(MAXX), B(MAXY), C(MAXX), D(MAXX) Assign dummy variables to the COMMON blocks NN = N DO I = I, N XX(I)= X(I) YY(1) = Y(I) END DO calculata H(I) DO 1=1 ,N- 1 H(r) = x(t+1)-X(I) ENDDO IF ( IBC EQ 1 ) THEN^Natural splines B(t)= I.D0 C(I) = O.DO D(1) = 0.D0 A(N) = ODO B(N) = I.D0 D(N) =0 DO ENDIF IF ( IBC .EQ. 2) THEN ! Clamped splines DERIVI = O.DO DERIV2 = 0.00 B(1) = 2D0 • H(I) C(1)= H(l) 13(1)= 3.D0*((Y(2)-Y(1))/H(1)-DERIVI ) A(N) = H(N- I ) B(N) = 2.D0•H(N-1) D(N) ---3.1)0*((Y(N)-Y(N- I))/H(N-1)-DERI VI ) ENDIF IF ( IBC .EQ. 3) THEN I Splines with fitted end points B(1)=-H(1) C(1)= H(I) D(1) = 3.D0*(H(I)**2) Y(1)/((X(1)-X(2)) • (X(I)-X(3)) • (X(I)-X(4))) + Y(2)/((X(2)-X(1)) • (X(2)-X(3)) • (X(2)-X(4))) + Y(3)/((X(3)-X(I)) (X(3)-X(2))* (X(3)-X(4))) Y(4)/((X(4)-X(I)) " (X(4)-X(2)) • (X(4)-X(3))) A(N)' H(N-I) B(N) --H(N-t) D(N)=-3.D0*(H(N-1)••2)* Y(N-3((X(N-3)-X(N-2))"(X(N-3)-X(N-I))*(X(N-3)-X(N))) ,Y(N-2)/((X(N-2)-X(N-3))"(X(N-2)-X(N-1))*(X(N-2)-X(N))) + Y(N-1((X(N-1)-X(N-3))"(X(N-1)-X(N-2))*(X(N-1)-X(N))) + Y(N)/((X(N)-X(N-3))*(X(N)-X(N-2))"(X(N)-X(N-I))) ENDIF A(1) DO DO I = 2, N-I IM =1-1 A(I) = H(IM) B(I) = 2 DO*(H(IM)+H(1)) C(I)= H(I) D(I) = 3.D0*((Y(1+1)-Y(1))/H(1)-(Y(1)-Y(IM))/H(IM)) END DO CALL TDMA ( A, B, C, D, R, I, N) DO I = I, N-1 IF =1+1 13(1) = (Y(1P)-Y(1))/H(1)-H(1)*(2 R(I)+R(IP))/3.D0 S(1) = (R(IP)-R(I))/(3.DO*H(I)) END DO RETURN END REAL'8 FUNCTION FS_INTEG ( Z I , Z2) IMPLICIT REAL*8 (A-H2O-Z) PARAMETER( MAXX=20I ) COMMON/ BLKA / X(MAXX), Y(MAXX), N COMMON/ BLKB / Q(MAXX), R(MAXX'), S(MAXX) IF (Z1.GT.Z2) PAUSE 'FEHLER' IF (ZI.LT.X(I)) THEN WRITE(", I 5) Z1 STOP 15^FORMAT(/'-Waming-',D10.5; is out of the range !!!!') ELSEIF (Z2.GT.X(N)) THEN WRITE(', 'I S) ZI STOP END1F C^Search for the start region of the Integral I=1 1=N 10 K=0+1)/2 IF(ZI.LT.X(K)) 1=K IF(Z1.GE X(K)) I=K 1F(1.GT.1+ I) GOTO 10 ISTART=1 C^Search for the end region of the integral 1=1 1=N 12^K=(I+1)/2 1F(Z2.LT.X(K)) 1=K 1F(Z2.GE.X(K)) I=K • 106 IF(1.GT.I+1) GOTO 12 C Calculate the integral over the functions FS_INTEG-0.D0 C^Now calculate the start K—ISTART KP=K+1 DX=X(KP)-X(K) FINT2=X(KP)*Y(K)+DX**2*(0.500*Q(K)+ • DX*0 .00/3.00*R(K)+ I . DO/4. DO* DX*S(K))) DX=ZI-X(K) FINTI =Z1* Y(K)+DX**2*(0.5DO*Q(K)+ • DX*0 .00/3.00*R(K)+ I .D0/4.00* DX* S(K))) FS_INTEG=FINT2-FINTI FS=1D0 DO K=ISTART+1,IEND KP=K+I DX=X(KP)-X(K) FINT2=X(KP)nY(K)+DX**2*(0.500.Q(K)+ DX*(1.00/3.D0*R(K)+1.00/4.DO*DX*S(K))) FINTI=X(K)*Y(K) FS INTEG=FS INTEG+FINT2-FINTI END 150 K=IEND KIP=K+ I DX=X(KP)-X(K) FINT2=X(KP)*Y(K)+DX**2*(0.5D0*Q(K)+ DX*(1.00/3.00*R(K)+1.00/4.00*DX*S(K))) DX=Z2-X(K) FINTI=Z2*Y(K)+DX**2*(0.500*Q(K)+ DX*0 .D0/3.DO*R(K)+1.00/4.DO*DX*S(K))) FS_INTEC,FS_INTEG-(FINT2-FINT1) RETURN END REAL*8 FUNCTION FS_INP ( Z ) C^This Function calculates the value for a spline interpolation. • The required values are supplied by a COMMON block_ C INPUT: • XARG: interpolation values for X • YARG: interpolation values for Y C^OUI-POT: • F2D: interpolation for Z(X,Y) IMPLICIT REAL•8 (A-H2O-Z) PARAMETER ( MAXX=201 ) COMMON! BLKA / X(MAXX), Y(MAXX), N COMMON! BLKB / Q(MAXX), R(MAXX), S(MAXX) IF(Z.LT.X( I)) THEN WRITE(*,15) Z 15^FORMAT(fl-Waming-',D10.5, is out of the range !!!!') ELSEIF (Z.GT.X(N)) THEN WRITE(*,-) Z I=N- I WRITE(*,15) Z ELSE 1=1 1=N 10^Kl+.1Y2 IFIZ LT.X(K)) 1=K IF(Z GE X(K)) I=K 1F(1.GT I+ 0 GOT° 10 ENDIF DX = Z - X(I) FS_INP=Y(1)+DX*(Q(IHDX*(R(1)+DX*S(I))) RETURN END HF2D_in.PUT Driver File ^ 101 ! ri of grids in X direction^[-) ^ II !8 of grids in R direction^[-] ^ 51 ! number of data in x direction in out-put file ^ 500.000 ! flow rate in the reactor^[ml/min] ^ 000.00 ! Inlet pressure^[Pa) ^ 00000.00 ! Outlet pressure^[Pa) ^ 21.500 - 2 ! fibre length ^ 115.000-6 ! Lumen radius^fro) ^ 124.000-6 ! Outer membrane radius^[ml ^ 175.00D-6 ! Krogh Radius^Em] ^ 0.0157500 inner shell diameter of the HFBR [m] ^ 8128 !6 of fibres in the reactor^I-) ^ 6.000-15 ! Membrane permeability Lp [m] ^ 5.000-15 ECS conductivity K^Dn21 ^ 10.00 ! BSA start-concentration^[gill ^ 6.8D-I0 ! BSA diffusivity^[rn2/s] ^ 69000.00 ! Molecular weight of the protein [g/moles] ^ 20.00 ! Temperature^toC/ ^ 0000 ! grid-spacing factor in x-dir. [-] ^ 0.000 ! grid-spacing factor in r-dir. [-] ^ 1.0D-6 ! Arc_ for the pressure^[Pa) ^ 1.900 ! SOR factor x direc. for L.b.L.(1 2) ^ 1.600 SOR factor in r direr, for L.b.L.(1 -2) ^ 1000 ^ ! max # of iteration to adjust pressures 1.0-8 ! steady-state acc for the concentrations .^! Time-step for the cliff cony. eq. [sl 2000 ^ 1.000 ! Time-relaxation factor for the duff, cony. eq. ^ 1.000 ! transient reporting values^[11) 2.0130 5.000 10.00 END 107 Source Code for the One-Dimensional Computer Model RETURN PROGRAM HFBR_Protein_lD VERSION 2.0 LAST CHANGE: 10.6.92 Numerical Solution for ID HFBR Polarization model! By setting the ECS conductivity to O. the empty ECS is simulated by: Jurgen Koska IMPLICIT REAL*8 (A-H2O-Z) CALL SET_UP START_TIME = F_T1ME_STAMP( ) CALL HEART CPU TIME = FTIMEETAMP( ) - START TIME WRITE(*,100) CPU_TIME/60. WRITE(20,100) CPU_TIME/60. CLOSE( UNIT=20 ) STOP 100 FORMAT( ' CPU-time: F12.4,' Iminj ') END SUBROUTINE HEART • Main Routine from where all other functions are called IMPLICIT REAL•8 (A-H2O-Z) PARAMETER ( MAXX=25I ) COMMON / B_DI / NX, NXM, RI, R2, R3, XL COMMON / N_OUTP / N_X_SKP COMMON / B_P I / P_LU(MAXX), P_OS(MAXX), PECS(MAXX) COMMON / C_I / C_RELAX, C_TIME_TOTA, C_DT, C_EPS COMMON / CONC / C_ECS(MAXX) COMMON/TROUT / NTRANS, TRATIME_OUT(99) COMMON / C_BEG1N / C_START COMMON / PLOT_TU IPLOT_TIME !PLOT = I N_FILE = 1 WRITE(20,125) 0, C_TIME_TOTA CALL PRESSURE ( ITER_P, 100, I .D-12 ) CALL DATA_OUT ( 0 , C_TIME_TOTA ) IF ( C_START .LE. 0.D0) RETURN CALL UP_DATE_C ( C_DT, TODO, C_DIFMAX, PE_X ) C TIME_TOTA = C_TIME_TOTA + C_DT DO ITER=I , 00000000 CALL UP_DATED ( C_DT, 0.5D0, C_DIFMAX, PE_X ) C_DIFMAX = C_DIFMAX / CDT CALL PRESSURE ( ITERP, 100, 1.D-12 ) C_TIME_TOTA = C_TIME_TOTA + C_DT TI_H = C_TIME_TOTA / 3600. ! convert to hours IF ( C_DIFMAX IE. CEPS ) THEN CC^the required steady-state criterion is obtained CALL DATA_OUT ( N FILE , TI_H WRITE(* ,200) WRITE(20,200) WR1TE(20,125) N_FILE, TI_H WRITE(20,130) ITER, C_DIFMAX, PEN RETURN END1F IF ( C_TIME_TOTA GE. TRA_TIME_OUT(N_FILE) ) THEN CC^We reached the time for writing all information WR1TE(20,125) N_F1LE, TI_H WRITE(20,130) ITER, C_D1FMAX, PE_X CALL DATA_OUT ( N_FILE , TI_H ) IF ( C_TIME_TOTA .GE. TRA_TIME_OUT(NTRANS) ) THEN CC^check if max. time is exceeded WRITE(20,300) WRITE(* ,300) RETURN END IF N_FILE = N_FILE + END IF IF ( C_T1ME_TOTA .GE. IPLOT_TIME "IPLOT ) THEN WRITE(*,135) ITER, TI_H, C_DIFMAX, PE_X CALL PLOT CALL DELAY ( LODO ) !PLOT = IPLOT + 1 END IF C DT = C_DT • C_RELAX EN-121 DO IOU FORM AT(5X,110,F I 6 ^ 2,1s1') 125 FORMAT(///' ^ file 6:15,' time. ,F9 4,' 130 FORMAT(2X,I5,6(2X,E12.4)) 135 FORMAT(15,2),F10.3,2X,F10.5,2X,F10.5,2X,F10.3) 200 FORMAT(//,10X,' steady state was reached !!!') 300 FORMAT(//,10X,' steady state was not reached !!!',//) END SUBROUTINE PRESSURE (ITER, MAXITER, EPS_P ) C Here the coupled lumen and ECS pressure eq. are solved for C convergence. IMPLICIT REAL'S (A-H2O-Z) PARAMETER( MAXX=251 ) COMMON / B_P I / P_LU(MAXX), P_OS(MAXX), P_ECS(MAXX) COMMON / B_004 / C_RELAX, C_T1ME_TOTA, C_DT, C_EPS COMMON / B_DI / NX, NXM, RI, R2, R3, XL CC first we update the osmotic pressures ! CALL UP_DATE_P_OS CC now we allow a max. number of iterations to calc. the pressures CC in the lumen and in the ECS for convergence. DO ITER = I, MAXITER CALL PRESS_ECS ( D_P_ECS ) CALL PRESS_LUMEN ( D_P_LU ) IF ( (D_P_LU+D_P_ECS) .LT. EPS_P ) THEN the pressures converged and we rale the velocities CALL VEL_ECS RETURN END IF END DO WRITE( *,100) C_T1ME_TOTA WR1TE(20,100) C_TIME_TOTA 100 FORMAT (2X,THE PRESSURES DID NOT CONVERGE AT:',F10.2,'s) CALL VEL_ECS PAUSE RETURN END SUBROUTINE PRESS_ECS ( P_D1FMAX ) C ^ now we calculate the new pressures in the ECS-side! IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=251 ) COMMON/ B_DI / NX, NXM, RI, R2, R3, XL COMMON / B_A / X JMAXX), D_X(0:MAXX) COMMON / B_SP / PER_MEM, COND_ECS COMMON / BPI / P_LU(MAXX), P_OS(MAXX), P_ECS(MAXX) COMMON / B_VISC / VISC_L, V1S_INS, V_K COMMON / C_1 / C_RELAX, C_TIME_TOTA, C_DT, CEPS COMMON / X BLK / X IN(MAXX), D_X_INV(MAXX) COMMON / PRESS / PIN , POUT COMMON / BLK_PS /-CO_PTECS,CO_P_LU, CO VEL DIMENSION A(MAXX), B(MAXX), C(MAXX), D(MAXX), P_NEW(MAXX) C^note: flux B.C. are given P_D1FMAX = 0.D0 A(I) = 0.D0 B(I) = - ( D_X_INV(I) + CO_P_ECS X_IN(1) ) C(I) =^D_X_INV(I) D(1) = ( - PIN - P_OS(I) ) • CO_P_ECS • XIN(l) DOI = 2, NXM IM = I - I A(I) = D_X_INV(IM) B(I) = -( D_X_INV(1M) + D_X_INV(1) + CO_P_ECS • X_IN(I)) C(I)=^D_X_INV(I) D(I) = ( - P_LU(I) - P_OS(I) ) • CO_P_ECS X_IN(I) END DO A(NX) = D_X_INV(NXM) B(NX)^D_X_INV(NXM) + CO_P_ECS • X_IN(NX) ) C(NX) = 0.D0 D(NX) = (- POUT - P_OS(NX) ) CO_P_ECS • X_IN(NX) CALL TDMA ( A, B, C, D, P_NEW, I, NX ) DO I = I, NX P_DIFMAX = DMAX I( P_DIFMAX, DABS(P_NEW(1)-P_ECS(I)) ) P_ECS(I) = P_NEW(I) 108 END DO RETURN END SUBROUTINE PRESS_LUMEN ( P_DIFMAX ) C^ This S.R. calculates the new lumen pressure profile! IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=251 ) COMMON / B_DI / NX, NXM, RI, R2, R3, XL COMMON / B_A / X_(MAXX), D_X(0:MAXX) COMMON / B_SP / PER_MEM , COND_ECS COMMON / B_PI / P_LU(MAXX), P_OS(MAXX), P_ECS(MAXX) COMMON / B_V1SC / V1SC_L, VIS _INS, V_K COMMON / X_BLK / X_IN(MAXX), D_X_INV(MAXX) COMMON / PRESS / P_IN, P_OUT COMMON / BLK_P_C/ CO_P_ECS,CO_P_LU, CO_VEL DIMENSION P_NEW(MAXX), A(MAXX), B(MAXX), C(MAXX), D(MAXX) P_DIFMAX = 0.00 C^ now we calculate the new axial pressures in the lumen with a C^ note: the pressures at both fibre ends are given as B.C. A(I) = 0.DO B(I)= 1.00 C(1) = O.DO D(I) = PIN ! pressure at the fibre lumen inlet DO I =2, NXM IM = I - I A(I) =^D_X_INV(IM) B(I) = - ( D_X_INV(1) + D_X_INV(IM) + CO_P_LU * X_IN(I) ) CO) = D_X_INV(I) D(I) = CO_P_LU • X_IN(I)* ( P_OS(I)-P_ECS(I) ) END DO A(NX)= 0.D0 B(NX)= 1.00 C(NX) 0.D0 D(NX)= P_OUT ! pressure at the fibre lumen outlet CALL TDMA ( A, B, C, D, P_NEW, I, NX ) DO I = I, NX P_DIFMAX = DMAX1( P_D1FMAX, DABS( P_NEW(4)-P_LU(1)) ) P_LU(I) = P_NEW(1) ! new pressure in axial direction END DO RETURN END • ...rt.*** ****** rt.*** ******** **.*****..*******.x* . ...... ......... SUBROUTINE SET_UP ........ ....*•* ..... ..... *******.*.”•.****• ............... CC This routine initializes the program by reading the external file CC and sets system parameters IMPLICIT REAL*8 (A-14,0-Z) PARAMETER ( MAXX=25I ) CHARACTER*8 FILE_NAME COMMON / F_NAME / FILENAME COMMON / B_PI / PI COMMON / B_A / X_(MAXX), D_X(0:MAXX) COMMON / B_PI / P_LU(MAXX), P_OS(MAX)C), P_ECS(MAXX) COMMON / C_I / C_RELAX, C_T1ME_TOTA, C_DT, C_EPS COMMON / B_DI / NX, NXM, RI, R2, R3, XL COMMON / B_SP / PER_MEM, COND_ECS COMMON/ HALF_GR / X_MID(MAXX+ I) COMMON / N_OUTP / N_X_SKP COMMON / B_VISC / VISC_L, VIS_INS, V_K COMMON / CONC / C_ECS(MAXX) COMMON / DIFFUS / DIFF, NFIBRE COMMON / OSMOTI / PM, SM, T, OS_CONS, VIR_I,VIR_2 COMMON / TROUT / NTRANS, TRATIME_OUT(99) COMMON / C_BEGIN / C_START COMMON / X_BLK / X_IN(MAX)G , D_X_INV(MAX)0 COMMON / BLK_PE / D_X_D1FF(MAXX) COMMON / PRESS / P_IN, P_OUT COMMON / BLK_P_C / CO_P_ECS, CO_P_LU, CO_VEL COMMON / PLOT_TI / 1PLOT_TIME OPEN( UNIT=10, FILE =1-1F1 D_IN.PUT) READ(I 0,333) FILE_NAME ! name for out-put files 333 FORMAT( A8 OPEN( UN1T=20, FILE = FILE_NAME/P.DAT) READ(I0,*) NX ! Of of grids in X-direction READ(10,") N_DAT 'skip factor for data in the out-put file READ(10,*) FLOW ! flow rate in the reactor [mL/min] READ(I 0,") PIN ! Inlet pressure [Pal READ(I0,") P_OUT 'Outlet pressure [Pal READ(I0,") XL^! [ml length of the capillary READ(I0,") RI^! [ml inner lumen radius READ(10,") R2^! [in] outer membrane radius READ(I0,") R3 ! lin] Krogh radius READ(I0,*) RSHELL '[ml inner shell radius of the HFBR READ(I0,") NFIBRE ! [-] (1 of fibres in the reactor READ(I0,*) PER_MEM ! lin) membrane permeability READ(10,*) COND_ECS! fin"2] ECS conductivity READ(10,•) C_START ! BSA start concentration in the ECS READ(I0,*) DIFF ! BSA diffusivity REA D(10,•1 PM^! molecular weight of the protein RE AD(10,*) TC^'temperature [KJ READ(10,*) G_FX^grid-spacing factor in x-dir. [-] READ(10,*) C_EPS ! steady-state acc. for the concentrations READ(I0,*) C _DT ! Time-step for the diff. cony. eq. REA D(10,*) IPLOT_TIME! Time-step for screen plots READ(10,") C_RELAX ! time-relaxation factor for the cliff, cone. eq. IF ( NX .GT. MAXX) NX=MAXX DO 1= I, 99 READ(I 0,",ERR=I8) TRA_TIME_OUT(I) END DO 18 N_TRANS = I - 1 CLOSE( UN1T=10 ) DOI= 1, N_TRANS TRATIME_OUT(I)= TRATIME_OUT(1)*3600.D0 ! convert to sec_ END DO VISC_L^= F_VISC_H20( TC ) V1S_INS^= 3.65D-3 V_K^= 2.085110 PI^= DATAN(I.D0) r 4.D0 = TC + 27115 C_TIME_TOTA = 0.110 NXM^= NX I N_X_SKP^= NXM / (N_DAT-1) CC^Check if the Krogh Radius is given)' IF ( R3 .LE. 0.00) THEN R3 = RSHELL / NFIBRE**0.5 WRITE() 'Krogh radius:', R3 END IF IF ( FLOW .GT. 0.000) THEN POUT = 0.000 QF = FLOW" I .D-6/60/NFIBRE 'calculate flow in a single fibre P_IN = QF' XL*8" VI SC_LJ(P1.121**4.) END IF WRITE(',100) P_1N, P_OUT 100 FORMAT(5)Qp-in:',F10.2,5X;p-out:',F I 0.2; [Pal') FLOW = P_IN/(I.D-6/60/NFIBRE" XL*8*VISC_L/(PI"R1*"4.)) WRITE(",110) FLOW 10 FORMAT(5X,'Q-HFBR:',F10.2,' [ml/minn CC-- Setup a non-uniform grid in both directions CALL USYM_GR1D (000, XL, NX, G_FX, X_, D_X ) DO I = I, NX P_LU(I) = FTLIN( 0.00, XL, P_1N, P_OUT, X(1) P_OS(I) = F_PI( C_START ) C_ECS(I) = C_START P_ECS(I) = (P_IN+P_OUT) "0.500 ! set start pressure X_MID(I) = X_(I) - D_X(1-1).0.5D0 END DO X_MID4NX+1) = X_(NX) P_LU(I) = P_1N P_LU(NX) = POUT DO 1= I, NX X_IN(I)^= ( D_X(I-1) + D_X(I) ) • 0.500 D_X_DIFF(I) = D_X(I ) / DIFF END DO DOI= I, NXM D_X_INV(I) = I.DO/D_X(I) END DO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC^check for the model and set the appropriate constants CC^K = 0 empty shell side model CC^K >0 filled ECS model with Darcy conductivity K CC^CO_P_ECS constant for pressures CC^CO_VEL constant for ECS velocities A_FRAC = 2. * RI / (R3**2.-R2**2.) IF ( COND_ECS .LE. O.DO ) THEN PHI = 4 *R3**4"DLOG(R3/R2)/(R3""2-R2"*2)-3"R3**2+R2**2 COND_ECS = PHI / 8 DO WRITE(',*) ' You use the empty ECS Model' END IF CO_VEL = - COND_ECS / V1SC_L CO_P_ECS = A_FRAC • PER_MEM / COND_ECS WRITE(*)' K= COND_ECS CALL DELAY ( 2.000 CC set constant for the lumen pressures CO_P_LU = 16.1)0 • PER_MEM / R 1" -3.110 CC setup constants for osmotic pressure relationship = 8.3141500 ! gas constant 109 VIR_ I = 0.01047300 ! 2nd canal coef for BSA VIR_2 = 1.7374313-5 ! 3rd virial coef for BSA OS_CONS = R • TI (PM/1000.) CALL DATA_SHOW( RSHELL, FLOW, GF_X ) RETURN END SUBROUTINE VEL_ECS Calculates the velocity profile in the ECS IMPLICIT REAL•8 (A-H2O-Z) PARAMETER ( MAXX=251 ) COMMON / B_DI / NX, NXM, RI, R2, R3, XL COMMON / B_A / X(MAXX), D_X(0:MAXX) COMMON! B_SP / PER_MEM, COND_ECS COMMON / B_P I / P_LU(MAXX), P_OS(MAXX), P_ECS(MAXX) COMMON! B_VE_ECS / VEL(MAXX,2) COMMON / B_V1SC / VISC_L, VIS INS, V_K COMMON / CONC / C_ECS(MAXX) COMMON! BLK_P_C / CO_P_ECS, CO_P_LU, CO_VEL CC^calculate velocities in x-direction U(x,r) e- and w-side VEL(I ,2)= 0.000 ! at X = 0 W DO I = 2, NXM, 2 IP =1+1 IM = I-I VEL(I,2) = CO_VEL " ( P_ECS(I)-P_ECS(IM) )/D_X(IM) VEL(1,1) = CO_VEL " ( P_ECS(IP)-P_ECS(I) )/D_X(1) VEL(IM,1) = VEL(1,2) VEL(IP,2)= VEL(1,1) END DO VEL(NX,1) = 0.0D0 ! at X = L E RETURN END SUBROUTINE PLOT •• Driver routine for graphical presentation on the screen IMPLICIT REAL*8 (A-HO-Z) PARAMETER( MAXX=251 ) EXTERNAL FS_1NP COMMON / B_D1 / NX, NXM, RI, R2, R3, XL COMMON! B_A / X(MAXX), D_X(0:MAXX) COMMON! CONC / C_ECS(MAXX) CA MAX =100,00 CALL SPLINE ( X_, C_ECS, NX, 2) CALL SCREEN ( FS_INP, 0.00, XL, 0.000, CA_MAX ) RETURN END B_N(I) = -( A_N(I) + C_N(I) + VEL(I, I ) - VEL(1,2) ) END DO A_N(N X) = PE_(NXM) + DMAX I (0.DO, VEL(NX,2) ) C_N(NX) =0.00 B_N(NX) = ( A_N(NX) - VEL(NX,2) ) DOI= 1, NX A(1) = F • A_N(1) C(1) = F • C_N(1) B(I) = F • B_N(I) - X_IN(I) R_DT 0(0= FM • A_OLD(I)^• C_ECS(I-1) + ( FM • B_OLD(I) - X_IN(1) • R_DT ) • C_ECS(I ) + FM • C_OLD(1)^• C_ECS(I+1) END DO CALL TDMA ( A, B, C, D, C_NEW, I, NX ) ! solve the tridiagonal set C save these values for the next time step DOI= I, NX A_OLD(I)= A_N(I) B_OLD(1) B_N(I) C_OLD(I)= C_N(I) C_D1FMAX = DMAX1( C_DIFMAX, DABS(C_NEW(1)-C_ECS(1))) C_ECS(1)= C_NEW(1) END DO RETURN END REAL•8 FUNCTION F_POWER ( PE) C^This function applies the power law scheme C^(Patankar, 1980) IMPLICIT REAL•8 (A-H, 0-Z) F_POWER= DMAX1( 0.DO, (1.00 - 0.100 * DABS(PE) )• '5.) END SUBROUTINE LUMEN_VEL IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=25I ) COMMON / B_PI / PI COMMON / B_D1 / NX, NXM, RI, R2, R.3, XL COMMON! B_A / X(MAXX), D_X(0:MAXX) COMMON! B_SP / PER_MEM , COND_ECS COMMON / B_VE_ME / VXR_LU(MAXX), UXR_LU(MAXX) COMMON! B_P I / P_LU(MAX)Q, P_OS(MAXX), P_ECS(MAXX) COMMON! B_VISC / V1SC_L, VIS_INS, V_K COMMON! DIFFUS /DIVE, NF1BRE DIMENSION VDUMMY(MAXX) C^ UXR_LU: contains the average axial velocities in the lumen C ^ VXR_LU: contains the radial lumen velocities at R2 C^ calculate the average lumen axial velocities SUBROUTINE UP_DATE_C ( DT, F, C_DIFMAX, PE_X_MAX ) CC This Routine solves for the new protein distribution in the ^ CC shell side of the HFBR. The transient diffusion convection eq. CC is solved with a modified Crank Nicolson Method. CC To ensure a tridiagonal dominant matrix the power law scheme CC with 1st order up-winding with is included. IMPLICIT REAL•8 (A-H2O-Z) PARAMETER( MAXX=251 ) COMMON! B_D1 / NX, NXM, RI, R2, R3, XL COMMON! B_A / X_(MAXX), D_X(0:MAXX) COMMON! B VE ECS / VEL(MAXX,2) COMMON! CON C- / C ECS(MAXX) COMMON! DIFFUS!137F F, NFIBRE COMMON! X_BLK / X_IN(MAXX), D_X_INV(MAXX) COMMON / BLK_PE / D_X_DIFF(MAXX) COMMON! BLK_OLD / A_OLD(MAXX), B_OLD(MAXX), C_OLD(MAXX) DIMENSION A_N(MAXX), B_N(MAXX), C_N(MAXX), PE(MAXX) DIMENSION A(MAXX), B(MAXX), C(MAXX), D(MAXX), C_NEW(MAXX) C_DIFMAX = 0.00 PE_X_MAX OD° FM^= - (I -F) R_DT = I .D0 / DT DO I = I, NXM PE^= VEL(I,I) • D_X_DIFF(1) PE (I) = F_POWER ( PE ) / D_X_DIFF(1) CC^calculate the highest local axial Pe number PE_X_MAX = DMAXI ( PE_X_MAX, PE) END DO CC note: VEL(i,l) velocity at east side of control volume CC^VEL(i,2) velocity at west side of control volume A_N(1) =000 C_N(I) = PE JD + DMAXI(0.D0,-VEL(1,11) B_N(1)^= (C_N(1)+ VEL(1,11) DO I = 2, NXM A_N(I)= PE(1-1) + DMAX1(0.D0, VEL(I,2) ) C_N(1) = PE(I) + DMAXI(O.D0,-VEL(1,1) ) DUMMY=. RI 2. 1(8." VISC_L) DO I =2, NXM 1P=I+ 1 IM = I. I SLOPE = (P_LU(IP)-P_LU(IM))/(D_X(1)+D_X(IM)) UXR_LU(I)= SLOPE DUMMY END DO UXR_LU(1) DUMMY • (P_LU(2)-P_LU(I) )/ D_X(I) UXR_LU(NX)= DUMMY * (P_LU(NX)-P_LU(NX-0)! D_X(NX-I) DUMMY = PER_MEM / V1SC_L DO I = I, NX VXR_LU(1) = DUMMY • ( P_LU(1) - P_ECS(I) + P_OS(I) ) VDUMMY(I) = DABS( VXR_LU(I) ) END DO Q_U_IN = UXR_LU(1) • PI " RI ••2. ! flow rate in the lumen C^calculate the leakage flow into the lumen CALL SPLINE ( X_, VDUMMY, NX, 2) O_V_LEK = 2.DOP1•121 • FS_INTEG ( X(I), X(NX) ) /2.00 FAC = I .D6"60."NFIBRE WRITE(20,55) Q_U_IN*FAC, O_V_LEK"FAC WRITE(• ,55) Q_U_IN•FAC,O_V_LEK•FAC RETURN 55^FORMAT(4X,'Q Lu.: ',F10.4; Q Leak: ',F10.4,' [mL/minf) END SUBROUTINE DATA_OUT ( NUM, TRANS:IT ME ) C ^ This Function writes the data to a external file IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=251 ) CHARACTER•3 1NT_CHAR, NUMBER CHARACTER•8 FILE_NAME COMMON / F_NAME / FILE_NAME COMMON! Bye_ME / VXR_LU(MAXX), UXR_LU(MAXX) COMMON! B_DI / NX, NXM, RI, R2, R3, XL COMMON! B_P1 / P_LU(MAXX), P_OS(MAXX), P_ECS(MAXX) ^ 110 COMMON / BA / XJMAXX) , D_X(0:MAXX) COMMON! CONC / C_ECS(MAXX) COMMON / B_VE_ECS / VEL(MAXX,2) COMMON / HALF_GR / X_MID(MAXX+1) COMMON / B_SP / PER MEM , COND_ECS COMMON / B_V1SC / VISC_L, VIS_INS, V_K COMMON / N_OUTP / N_X_SKP 101^= I I NUMBER = INT_CHAR( NUM) OPEN( UNIT=I01, FILE = FILE_NAME/MNUMBER DO 1 = I, NX IF ( C_ECS(I) .LT. 0.00 ) C_ECS(I)= 0.00 END DO CALL LUMEN_VEL ! calculate the lumen velocities CALL C_AVERAGE ! calculate a mass-balance over the reactor DO I = 1, NX, N_X_SKP IM = I - I WRITE(101,100) X_(I) 1.D2, C_ECS(I), VXR_LU(1), interpolate the axial velocity in the ECS at the grid-points (D_X(1)•VEL(1,2)+D_X(IM)*VEL(1,1))/(D_X(1)+D_X(IM)), UXR LU(I), PLUM, P_ECS(1),P_OS(1),PECS(1)-P_OS(1) END DO WRITE(101,110)TRANS_TIME! time-stamp on each file CLOSE( UNIT=I01 ) RETURN 100 FORMAT( 2(F10.4,2X),3(E16.8,2)0,4(F10.4,2X) ) 110 FORMAT ( /F14.6) END SUBROUTINE C AVERAGE ............ .. 11■•••• .• •.* .................^..... ••••■■•• This S. R. calculates the average concentration in the HFBR IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=251, MAXR=I I) COMMON / B_A / X JMAXX), D_X(0:MAXX) COMMON! B DI / NX, NXM, RI, R2, R3, XL C / C ECS(MAXX) COMMON / CON COMMON / C_BEG1N7 C_START CALL SPLINE ( X_ ,C_ECS, NX, 2) CAV = FS_INTEG ( X(I), X JNX))/ X (MO WRITE(20,100) CAV , C_START-CAV WRITE( •,100) CAV , C_START-CAV 100 FORMAT(5X,'mass-balance:',F8.4,' lost.',F10.6,' [g/Lr) RETURN END X.* ..... ••■■•*. ............ 41.****.**.*. .....^ ..... rtt• • .***• SUBROUTINE DATA_SHOW ( RSHELL, FLOW, GF_X ) •••• ********** ** ****************** ••• •• ••• •• • ********* "`• • •• • •• •• All system information is written loan external file! IMPLICIT REAL•8 (A-H2O-Z) PARAMETER ( MAXX=251 ) CHARM. I ER*8 FILE NAME COMMON / F_NAME-/ FILE NAME COMMON! B_DI / NX, NXM, RI, R2, R3, XL COMMON / N_OUTP / N X SKP COMMON / B_A /^D X(01MAXX) COMMON / C_I / C RELAX, C TiME _TOTA, CDT, CEPS COMMON / B SP / PER_MEM, COND_ECS OMMON / DIFFUS 1fF, / NFIBRE D C COMMON / OSMOTI / PM, SM, T, OS CONS, VIR 1,VIR_2 COMMON / TROUT / N TRANS, TRATIMEDUT(99) COMMON / C^/ J_START COMMON! PRESS / P_IN, P_OUT COMMON / B_PI / PI 10 = 33 OPEN( UNIT=I0, FILE = FILE_NAME/P.INF) WRITE(I0,001) NX 001 FORMAT(6X;x-grids: ',I4) WRITE(10,003) P IN •I.D-5 ,P_OUT•I.D-5 003 FORMAT(2X,'BMss ures- inlet^:T10.4,' [bur outlet^[burl') WRITE(I0,005) FLOW 005 FORMAT(2X;Vol flow rate in HFBR ^.',F10.4,' [mL/minr) WRITE(10,007) XL•100. 007 FORMAT(2X,'Fibre length^:',F10.4,' [cmj') WRITE(I0,009) RI•106,R2` I D6 009 FORMAT(2X,'Fibre radii- inner^:',F10_4,' 4u-nr/ outer^:',F10.4,' [pm[') WRITE(I0,011) R3•1D6 011 FORMAT(2X,'Krogh radius^:',F10.4; lima) WRITE(I0,013) RSHELL • 100. 013 FORMAT(2X,'Inner cartridge radius^:',F10.4, [cm]) WRITE(I0,015) NFIBRE 015 FORMAT(2X;# of fibres in reactor ^:',I10) WRITE(I0,016) P1121•2XLNFIBRE Pl•R2"2•XL*NFIBRE 016 FORMAT(2Vinner HFBR membrane area :',F10.4; [inn/ 2X;outer HFBR membrane area^:',F10.4: [m2[') WRITE(I0,017) F_VISC_H20( T - 273.15) 017 FORMAT(2X,'Viscosity in lumen ^:',E10.4; [Pa-sr) WRITE(I0,019) PER_MEM 019 FORMAT(2X,'Membrane permeability ^:',E10.4,' [ml') WRITE(I0,020) PER_MEM • RI * DLOG (R2/R1) 020 FORMATRX,'calc. Membrane conductivity :',E10.4; [m2r) WRFTE(10,021) CON D_ECS 021 FORMAT(2X,'ECS-conduchvity ^:',E10.4,' frn2r) WRITE(10,023) 023 FORMAT(/,10X; Protein:) WRITE(I0,025) C_START 025 FORMAT(2X,'Stan Concentration ^:',E10.4,' [g/L[') WRITE(10,027) DIFF 027 FORMAT(2X,'Diffiisivity^:',E10.4,' [m/s••21) WRITE(I0,029) PM 029 FORMAT(2X,'Molecular weight^[g/molesr) WRITE(I0,031) F_PI(C_START) 031 FORMAT(2X,'Osm. press. ^:',E10.4,' [Pal') WRITE(10,032)F_VISC_BSA(C_STARD 032 FORMAT(2X;Viscosity of protein sol. :',E10.4,' [Pa-sr) WRITE(10,033) VIR_I ,VIR_2 033 FORMAT(2X,'I at virial coeff. ^:',F12.8/ .^,2X,'2 nd vinal coett^:',E12.4) WRTTE(I0,041) T-273.1500 041 FORMAT(2X,Temperature^:',F I 0.2; [oCr) WRITE(I0,043) C_EPS 043 FORMAT(2X,'Steady-state for conc.^:',E10.4,' [g/L,]') WRITE(10,045) C_DT,C_RELAX 045 FORMAT(2X,Time Step ^:',F10.4,/ 2 X,Tirne Aced. Factor:^:',F10.6) WRITE(10,049) 049 FORMAT(//,2X,'"• • Transient times for mint-out ••• ') DOI= I, NTRANS TR = TRA_TIME_OUT(I) WRITE(I0,055) TR, TR/60.,TR/60./60. END DO 055 FORMAT(5X, P102,' [sj^[min] ',F10.4,' [h]) WRITE(I0,057) GF_X 057 FORMAT(//,2X,'•*"' grid locations in x-ditection [cm] "'"• ',F10.6) DO I = 1, NX WRITE(I0,059) 1,X J1)*100., D_X(05100. END DO 059 FORMAT(2X,I4,2(5X, P10.6)) 063 FORMAT(2X,I42(5X, P10.6)) CLOSE( UNIT=I0 ) RETURN END ***** *****^********* ..***.rt ******^****** *4..**.***. SUBROUTINE UP_DATE_P_OS CC This Sub. R. updates the osmotic pressures along the membrane/ECS CC interface. IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=251 ) COMMON! B_DI / NX, NXM, RI, R2, R3, XL COMMON / BPI / P_LU(MAXX), P_OS(MAXX), P_ECS(MAXX) COMMON / CONC / C_ECS(MAXX) DO I = 1, NX P_OS(I)= F_PI( C_ECS(I) ) END DO RETURN END REAL•8 FUNCTION F_VISC_H20( TC ) C This function calculates the viscosity of H2O at the given Temp TC in Celsius. C input IC^oC C output vise^Pa s IMPLICIT REAL•8(A-H,0-Z) F_VISC H20 = 10.D0• *( (1.3272.(20-TC)-0.001053*(TC-20)""2 ) 7(TC+105.) ) • 1.0021/100/10 RETURN END REAL*8 FUNCTION F_VISC_BSA( C) C This function relates the viscosity to the BSA concentration. C The viscosity of scaler must be defined before in VISC_L C input c BSA g/L C output vise Pa s IMPLICIT REAL•8(A-H2O-Z) PARAMETER ( MAXX=251 ) COMMON / B_ VISC / VISC_L, VIS_INS, V_K C2 = C/10 F_VISC_BSA = VISC_L•( 1 • C2 * ( VIS INS + V_K • C2 ) ) RETURN END 111 REAL*8 FUNCTION F_PI( C_IN ) ^ CC This function calculates the local colloid osmotic pressure as a ^ CC function of the local protein concentration according. ^ IMPLICIT REAL*8 (A-H2O-Z)^ COMMON / OSMOT I / PM, SM, T, OS_CONS, VIR_I,VIR_2 ^ F_PI OS_CONS C_IN •^ . (I + C_IN • ( VIR_1 + C_IN • V1R_2 ) ) ^ RETURN^ END^ REAL*8 FUNCTION F_YL1N( XI, X2, Y I , Y2, XLIN )^ this Function performs a linear interpolation between two points. ^ IMPLICIT REAL*8 (A-H2O-Z)^ F_YLIN ( Y2 - YI )/( X2 - XI)' ( XLIN XI) + Y1 RETURN END^ CHARACTER*3 FUNCTION 1NT_CHAR (I) ^ C This Function changes a an integer loan char. string (0-99)^ IMPLICIT REAL*8 (A-H2O-Z)^ CHARACTER*1 Al, A2, A3^ ITEL = I/10^ Al ='0'^ A2 = CHAR( 1TEL + 48 )^ A3 = CHAR( I - rret..lo 48 )^ INT_CHAR = A I //A2//A3^ RETURN^ END^ SUBROUTINE USYM_GRID ( XO, XF, N, G_F, X, DX)^ C This S.R. establishes a non uniform grid.^ C For positive G_F values more grids are on the left side. ^ C For negative G_F values more grids are on the right side. ^ C^0 < GF_A < PI/2 IMPLICIT REAL*8 (A-H2O-Z) DIMENSION X(N), DX(0:N)^ PIH = 2.D0 * DATAN(I .D0) IF ( DABS(G_F) .(JT. PIH)^ NM = N - I^ IF (0_F .EQ. 0.00) THEN^ DXC = (XF-X0)/NM X(t)= X0^ DO I = 2, N^ X(I) = X(1-0+ DXC^ DX(11- I DXC END DO ELSE^ G = DABS(G_F)^ DO I = I, N^ 1M — I - I^ X(I) = XO + (XF-X0) • DIAN( 0 IWNM )/DTAN(G)^ IF ( I .NE. I) DX(IM )= X(I)-X(IM)^ END DO IF (0_F .LT. 0.00) THEN^ DO I = I, N^ DX(I) = X(N-I+1) - X(N-I)^ END DO^ X(I) = X0^ DO I = 2, N^ X(I) = X(1-0+ DX(I-1)^ END DO^ END IF^ END IF^ DX(0) = 0.00^ DX(N) = 0.D0^ RETURN^ END^ SUBROUTINE SYM_GRID ( XO, XF, N, 0_F, X, DX)^ IMPLICIT REAL*8 (A-I-1,0-Z) C This S.R. establishes a symmetric non uniform grid C^0 < GF_A < P1/2 C^Note GF_A = 0.00 leads to a uniform grid DIMENSION X(N), DX(ON) PIH = 2.00 * DATAN(I .D0) IF ( DABS(G_F) .GT PIH) NM = N - 1 IF ( G_F .EQ. 0.D0 ) THEN DXC = (XF-X0)/NNI X(1)= X0 DO I = 2, N X(I) = X(I-1) + DXC DX(1-1)=- DXC END DO ELSE 0 = DABS(G_F) DO I = I, N 1M -- I - I X(I) X0 + (XF-X0) • 0.5D0 " (I .DO+DSIN( G " (2.D0 • IM/NM-1.D0 ))/DSIN(G) ) IF ( I NE. I ) DX(IMX(1)-X(IM) END DO END IF DX(0) = 0.D0 DX(N)= O.DO RETURN END SUBROUTINE TDMA ( A, B, C, D, X, NI, N2) CC This S.R. solves a tridiagonal system of equations. CC In this routine the start and end of the eq. in the CC matrices has to be defined. IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=25I ) DIMENSION A(N2), B(N2), C(N2), D(N2), X(N2), P(MA30(), Q(MAX)) P(NI) = -C(NI)/ B(N1) COD = D(NI)/B(NI) DOI = NI+1,N2 IM =1-1 DEN = A(I)* POM)+ B(I) P(I) = -Ca) / DEN Q(I)=- (1)(0-A(1)•Q(IM))/DEN END DO X(N2) = Q(N2) DOI =N2-1,N1, -I X(I) = P(I)* X(I+1) + Q(I) END DO RETURN END SUBROUTINE SPLINE ( X, Y, N, IBC) C^This subroutine calculates the Q,R,S values for the spline interC^polation. C^The boundary for each end can be selected out of the three types: C^I: natur-al splines C^2: clamped splines C^31 fitted splines IMPLICIT REAL*8(A-H2O-Z) PARAMETER ( MAXX=251 ) COMMON / BLKA / XX(MAXX),YY(MAXX), NN COMMON / BLKB / Q(MAXX), R(MAXX),S(MAXX) DIMENSION X(N), Y(N), H(MAXX) DIMENSION A(MAXX), B(MAXX), C(MAXX), D(MAXX) C Assign dummy variables to the COMMON blocks NN = N DO I = I, N XX(I) = X(1) YY(I)= Y(I) END DO C^calculate H(I) DO I=1,N-1 H(0= X(I+ 0-X(1) ENDDO IF ( IBC .EQ. I ) THEN !^Natural splines B(I)= I.00 C(I) = 0.00 D(1) = 0.D0 A(N)= O.DO B(N) = I .D0 D(N) = 0.00 ENDIF IF ( IBC .EQ. 2) THEN ! Clamped splines DERIVI =000 DERIV2 -0.00 B(1) = 2.D0 • H(1) CID= H(I) D(I) = 3.D0*((Y(2)-Y(1))/H(1)-DERI VI ) A(N) = B(N) = 2 DO•1-1(N -I) D(N)=-3.00"((Y(N).Y(N- I ))JH(N-1)-DER1 V 1 112 ENDIF IF ( IBC .EQ. 3 ) THEN ! Splines with fitted end points B(I ) -H(I) C(I)= H(I) D(I) = 3. D0-(H(1)**2) • Y(1)/0X(1)- X(2)) " (X(1)-X(3)) (X(1)-X(4))) Y(2)/((X(2)-X(I)) (X(2)-X(3)) • (X(2)-X(4))) + Y(3)/((X(3)-X(1)) • (X(1)-X(2)) • (X(3)-X(4))) Y(4)/((X(4)-X(I)) • (X(4)-X(2)) • (X(4)-X(3))) A(N) = H(N-1) 13(1,1) =-H(N-1) D(N)=-3.D0*(H(N-1)•"2)" Y(N-3)/((X(N-3)-X(N-2))*(X(N-3)-X(N-1))*(X(N-3)-X(N))) + Y(N-2)/((X(N-2)-X(N-3))*(X(N-2)-X(N-1))*(X(N-2)- X(N))) + Y(N- I y((c(N- D-xos,-3))*(x04-0-x(N-2))*(x(N-1)-X(N))) Y(N)/((X(N)-X(N-3))*(X(N)-X(N-2))•(X(N)-X(N-1))) ENDIF A(I DOI= 2, N- I IM =1-I A(I)= H(IM) B(I) = 2.D0*(H(IM)+H(1)) C(I)= H(1) D(I) = 3 .D0•((Y(1+ 1)- Y(I))/H(I)-(Y(I)- Y(IM))/H(IM)) END DO CALL TDMA ( A, B, C, D, R, 1, N) DO! = I, N-1 IP =1+1 Q(I) = (Y(IP)-Y(1))/H(1)-H(1)*(2.*R(1)+R(IP))/3.D0 S(I) = (R(IP)-R(I))/(3.DO•H(I)) END DO RETURN END K=IEND KP=K+I DX=X(KP)-X(K) FINT2=X(KP)*Y(K)+DX**2*(0.500•Q(K)-+ DX*0 .D0/3.DO"R(K)+1.00/4.DO•DX•S(K))) DX-'Z2-X(K) FINT I =Z2" Y(K)+DX**2*(0.5DO*Q(K)+ DX*(1.1,10/3.DO*R(K)+1.D0/4.DO•DX•S(K))) FS_INTEG=FS_INTEG-(FINT2-FINTI) RETURN END ............. * ...... * ....... * .......... ***********•** ****** ***** REAL*8 FUNCTION FS_INP ( Z ) This Function calculates the value for a spline interpolation. The required values are supplied by a COMMON block. INPUF: XARG: interpolation values for X YARG: interpolation values for Y OUTPUT: F2D: interpolation for Z(X,Y) IMPLICIT REAL*8 (A-H2O-Z) PARAMETER ( MAXX=251 ) COMMON / BLKA / X(MAX)Q, Y(MAXX), N COMMON / BLKB / Q(MAXX), R(MAXX), S(MAKX) I=I J=N ^ 10 K = (1+4)/2 IF ( Z .LT. X(K) ) J=K IF ( Z .GE. X(K) ) I=K IF ( J .GT. 1+1 ) GOTO 10 DX = Z - X(I) FS_INP=Y(1)+DX*(Q(1)+DX*(R(1)+DX•S(1))) RETURN END REAL•8 FUNCTION FS_INTEG ( Z1, Z2) IMPLICIT REAL•8 (A-H2O-Z) PARAMETER( MAXX=251 ) COMMON/ BLKA / X(MAXX), Y(MAXX), N COMMON/ BLKB / Q(MAX50, R(MAXX), S(MAXX) IF (Z1.GT.Z2) PAUSE 'FEHLER' IF (ZI .LT.X(1)) THEN WRITE(`,15) Z I PAUSE 'FEHLER' 15^FORMAT(5X,'Waming :,E10.5; is out of the range !!!!') ELSE1F (Z2.GT.X(N)) THEN WRITE(•,15) ZI STOP ENDIF C^Search for the start region of the integral 1=1 .1=N 10^KI+J)/2 IF(Z1.LT.X(K)) J=K IF(ZI.GE.X(K)) I=K IF(J.GT.1+1) GOTO 10 ISTART=I C^Search for the end region of the integral 1=1 J=N 12^Kl+.4)/2 IF(Z2.LT.X(K)) J=K IF(Z2.GE.X(K))1=K IF(.1.GT.1+1) GOTO 12 IEND = I C^Calculate the integral over the functions FS_INTEG.D0 C^Now calculate the start K=ISTART KP—K+1 DX=X(KP)-X(K) FINT2=X(KP)*Y(K)+DX••2*(0.5D0•Q(K), DX*(1.00/3.DO•R(K)+1.D0/4.DO•DX"S(K))) DX=Z1-X(K) FINTI=ZI*Y(K)+DX••2*(0.5D0•Q(K)+ DX•O.D0/3.13.0*R(K)+1.D0/4.DO•DX•S(K))) FS_INTEG=FINT2-FINT1 FS0.D0 DO K=ISTART+1,IEND KP=K+1 DX=X(KP)-X(K) FINT2=X(KP)Y(K)+DX••2•(0.5DO•Q(K)+ DX*(1.D0/3.DO•R(K)+1.D0/4.DO•DX•S(K))) FINT I =X(K)*Y(K) FS_INTEG=FS_INTEC+ F INT2-F INT I END DO SUBROLTTINE DELAY ( TIME_STEP ) .* .....^ .....^..... ....... IMPLICIT REAL•13(A-H2O-Z) START_T = F_TIME0 100 IF ( ( F_TIME0 START_T) .GT TIME_STEP ) THEN RETURN ELSE GOTO 100 END IF END REAL*8 FUNCTION F_TIME_STAMPO C This function returns the time in seconds after midnight. IMPLICIT REAL*8 (A-H2O-Z) INTEGER"2 hour, minute, second, hundredth, YEAR, MONTH, DAY CALL GETT1M( hour, minute, second, hundredth ) CALL GETDAT (YEAR, MONTH, DAY) WRITE(,100) DAY,MONTH,YEAR, HOUR, MINUTE, SECOND F_TIME_STAMP = ((DBLE( hour )"3600.)+(DBLE( minute)*601+ DBLE( second) + (DBLE( hundredth ) / 100.0)) 100 FORMAT (/,2X:date: .^/,2X:time: RETURN END REAL*8 FUNCTION F_TIME( ) C This function returns the time in seconds after midnight. 1NTEGER•2 hour, minute, second, hundredth CALL GETTIM( hour, minute, second, hundredth ) F_T1ME = ((DBLE( hour ) • 3600.) + (DBLE( minute) • 60)' .^DBLE( second) + (DBLE( hundredth ) / 100)) RETURN END SUBROUTINE SCREEN ( FUNC, XI, X2, YSIAL, YBIG ) CC Performs a graphical screen print IMPLICIT REAL•8 (A-H2O-Z) PARAMETER ( ISC R=65 ,JSCR=22 ) CHARACTER"' SCR(ISCR,JSCR), BLANK, ZERO, YY, XX, FF DIMENSION Y(ISCR) DATA BLANK, ZERO, YY, XX, FF/ DO 1= 1, JSCR SC R(1,^= YY SCR(ISCR,J)= YY END DO 113 DO I = 2, ISCR-I SCR(l, I) = XX SCR(1, JSCR) = XX DOS = 2, 50CR-1 SCR(I,J)= BLANK END DO END DO DX = (X2 - XI )! (ISCR - I) X = XI YBIG = 0.D0 C YSML = YBIG DO! = 1, ISCR Y(I) = FUNC(X) C^IF ( Y(I) .LT. YSML ) YSML = Y(I) IF ( Y(I) .GT. YBIG ) YBIG = Y(I) X=DX+X END DO IF (YBIG .EQ. YSML) YBIG=YSML+ I DYJ = (JSCR-I)/ (YBIG-YSML) JZ = 1 - YSML • DYJ DO! = I, ISCR SCR(1,1Z) = ZERO = 1 + ( Y(1)- YSML ) * DYJ SCR(1,1)= FF END DO WRITE(,100) YBIG/10,(SCR(1,1SCR),1=1,1SCR) DO = JSCR-I, 2, -I WRITE(•,110)(SCR(1,1),1=1,1SCR) END DO WRITE(*,100) YSMU10,(SCR(1,1),1=1,ISCR) WRITE(*,I20) X1/10, X2110 100 FORMAT(IX,IPF10.3,1X,80A1) 110 FORMAT(1 2X,80A1) 120 FORMAT(' OX, I PF63,55X,F7.3) RETURN END Driver file "HFIDJN.PUT" REQUIRED the ID HFBR Protein Polarization Program ^ Test ! define name for output files ^ ^ 201 of grids in X direction I-1 ^ 101 O of data in the out-put file ^ 500.0D0 ! flow rate in the reactor 1mL/nunl ^ 0000.D0 ! Inlet pressure Raj ^ 000.0D0 ! Outlet pressure [Pal ^ 0.215D0 ! Length of the capillary ^Ern 115.00D-6 ! Lumen radius [ni ^ 124.00D-6 ! Outer membrane radius ^fm 175.00D-6 ! Krogh Radius fm 0.01575D0 ! inner shell diameter of the HFBR Ins ^ 8128 of fibres in the recator ^ 6.00D-I 5 ! Membrane permeability Lp ^Em 5.00D-15 ! ECS conductivity K [m••21 ^ 10.0D0 ! BSA start-concentration [0] ^ 6.8D-I0^! BSA diffusivity fin2/si 69000.D0 ^! Molecular weight of the protein Eg/molesj 20.DO ^! Temperature^foci 0.0D0 ^ ! grid-spacing factor in x-dir. ^f-1 I D-8 ^! steady-state acc. for the concentrations 1000.D0 ^! Time-step^ fs1 50000 ^ !Time-step for screen plots^[s] 1.000D0 ! Time-relaxation factor for the cliff cony. eq. 1000.D0 ! time request for transients 2000.D0 END
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Protein polarization in packed hollow fibre bioreactors
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Protein polarization in packed hollow fibre bioreactors Koska, Jurgen 1993
pdf
Page Metadata
Item Metadata
Title | Protein polarization in packed hollow fibre bioreactors |
Creator |
Koska, Jurgen |
Date Issued | 1993 |
Description | Osmotically active proteins and cells, retained in the extra capillary space (ECS) of hollow fibre bioreactors (HFBRs), can influence the hydrodynamics of such devices. A mathematical model was developed to describe the coupled hydrodynamics and high molecular weight protein transport in a cell filled HFBR. It was assumed that the multi-fibre reactor can be represented by a single, straight fibre surrounded by a symmetry envelope containing fluid and a homogeneous packed bed of cells. The low Reynolds number flow in this porous medium was described by Darcy's law. Because of difficulties associated with operating a reactor filled with mammalian cells, a suitable analogue was used to experimentally investigate protein polarization in packed HFBRs. The ECS side of the reactor was filled with an agarose/protein solution which, upon cooling, formed a porous medium with a uniform initial concentration of protein. A constant lumen flow was established for several days before the cartridge was sacrificed and the axial ECS protein distribution was measured. Since the protein transport and HFBR hydrodynamics were coupled, numerical methods were required to solve the governing equations of both the two-dimensional (axial plus radial variations) and the one-dimensional (axial variations only) models developed to predict axisymmetric transient ECS protein concentrations. Computer modelling results indicated that, because of the large length/radius ratio of the representative fibre unit, the two-dimensional ECS protein concentrations could be accurately duplicated by the simpler one-dimensional model. The latter model, required about two orders-of-magnitude less computational time than the former. The one-dimensional model results were compared to experimentally obtained ECS protein profiles and, subsequently, the model was used to predict protein polarization in ITFBRs for different conditions. The hydraulic conductivity of the agarose gel, required for the model, was experimentally determined using the falling head method. The measured conductivity values failed to adequately describe the observed protein polarization in the ECS of HFBRs. However, by using a gel conductivity which was about an order-of-magnitude higher than the measured value, it was found that the model agreed well, in general, with the measured ECS protein polarization profiles obtained for initial protein loadings of 5 - 30 g/L. The higher apparent conductivity, needed to fit the model to the experimental results, was attributed to the inability of the gel to completely fill the geometrically complex ECS. Since in HFBRs, packed anchorage-dependent mammalian cells are expected to achieve hydraulic conductivities similar to values encountered in tissues, several model simulations were carried out at low tissue conductivity values. The results indicated that, for these conditions, protein transport in the ECS is mainly governed by diffusion. Protein polarization, a dominant feature of empty ECS protein transport, is greatly reduced. Also, it was shown that the removal of product protein from a packed ECS space can be difficult, since ECS flows are reduced. Models, such as those developed here, can be used to further investigate HFBR operation and process optimization. |
Extent | 5278933 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2008-09-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0058557 |
URI | http://hdl.handle.net/2429/2043 |
Degree |
Master of Applied Science - MASc |
Program |
Chemical and Biological Engineering |
Affiliation |
Applied Science, Faculty of Chemical and Biological Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1993-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1993_fall_koska_jurgen.pdf [ 5.03MB ]
- Metadata
- JSON: 831-1.0058557.json
- JSON-LD: 831-1.0058557-ld.json
- RDF/XML (Pretty): 831-1.0058557-rdf.xml
- RDF/JSON: 831-1.0058557-rdf.json
- Turtle: 831-1.0058557-turtle.txt
- N-Triples: 831-1.0058557-rdf-ntriples.txt
- Original Record: 831-1.0058557-source.json
- Full Text
- 831-1.0058557-fulltext.txt
- Citation
- 831-1.0058557.ris
Full Text
Cite
Citation Scheme:
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>
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.831.1-0058557/manifest