Mathematical modeling of the batch adsorption of proteins on new restricted access media with poly(ethylene glycol) as a semi-permeable barrier using compact finite differences

Abstract

In this work we describe a mathematical analysis of the batch adsorption process of several proteins using a new restricted access medium consisting of agarose beads grafted with poly(ethylene glycol) (PEG) as a semi-permeable barrier and immobilized metal ions or ion exchange groups as binding sites. The model was fitted to experimental data, allowing the estimation of the adsorption rate constant and the effective diffusivity for each protein. The model was solved using compact finite differences in a MATLAB® platform. According to the results, the presence of grafted PEG reduces the adsorption of all proteins to different extent; with high molecular weight proteins being affected the most. The model also establishes a reduction in the adsorption rate constant (which affects protein interaction with binding sites). The movement of the protein molecules in the adsorbent pores is also affected by the grafted PEG, but to a lesser extent.

Keywords:
Restricted access medium; Batch adsorption modeling; IMAC and IEX; Compact finite differences

INTRODUCTION

The mathematical description of the batch protein adsorption phenomenon is always required since it allows estimating the adsorption rate constant, the effective diffusivity, the maximum protein capacity, and the desorption equilibrium constant for a particular protein with a particular adsorbent. These parameters can then be used to predict the behavior of protein adsorption using fixed beds.

The maximum protein adsorption capacity and desorption equilibrium constant can be estimated from adsorption isotherms (equilibrium data), while the adsorption rate constant and the effective diffusivity must be determined from kinetic adsorption studies. Thus, these parameters not only allow describing and elucidating the interaction between the target protein and the adsorption sites, but can also help characterize column chromatographic experiments (Sharma and Agarwal, 2002; Gutierrez et al., 2007).

In the present analysis, a generic three-resistance model where the adsorption of proteins from the bulk liquid to the adsorbent solid involves resistances to mass transfer such as film diffusion, pore diffusion, and surface reaction kinetics was implemented to describe the batch adsorption uptake of proteins. The model is based on basic physical and mass balance principles along with appropriate constitutive equations (e.g., second order forward reaction kinetics) and has been previously implemented to describe the adsorption, in a finite bath, of β-galactosidase onto anti- β-galactosidase immobilized on porous silica (Arve and Liapis, 1987). The same three-resistance model was also used to describe the affinity adsorption of immunoglobulin G to protein A immobilized on porous agarose matrices (Horstmann and Chase, 1989). The model has also been used to characterize the adsorption process in a column packed with spherical adsorbents (Carlsson et al., 1994; Kempe, 2007). Simplified models have been used to describe the adsorption of proteins onto spherical adsorbents. Skidmore et al. (1990) applied a lumped parameter model in which only the surface reaction rate resistance is considered relevant to describe the adsorption of protein using affinity adsorbents. In the same way, Arnold et al. (1985) presented a model that considers all the mass transfer resistance associated only to pore diffusion. These last two models, as a result of their simplicity, permit one to obtain an analytical solution of the batch adsorption process; however, not quite general.

Since the three-resistance model described in this work is composed of non-linear coupled partial differential equations, an analytical solution is not possible and, therefore, an implementation of a numerical solution is required. To accomplish this task, two approaches can be considered, the first method consists in the discretization of the time and space radial coordinates to generate a system of equations that needs to be solved simultaneously in each time step. The second one, called the method of lines and used in the present work, discretizes only the radial coordinate to generate a system of ordinary differential equations integrated in the time coordinate using a Runge-Kutta method or any other equivalent numerical approximation. For discretization in the radial coordinate, in the present work, seldom used compact finite differences (Lele, 1992) were employed. Using a three-point stencil, a fourth-order scheme can be developed with less truncation error when compared to conventional finite differences (Gurarslan, 2010). The parameter estimation, of unknown model parameters, can also be easily performed using a MATLAB® routine.

The adsorbents evaluated in this work were restricted access chromatographic media which consist of agarose-based adsorbents modified with adsorption sites (ion exchange groups or immobilized metal ions) and with grafted PEG as the controlled permeation agent. These adsorbents are able to adsorb small molecules (peptides or proteins) under the presence or large proteins which are prevented from reaching the adsorption sites (Gonzalez-Ortega et al., 2012a,b).

THEORETICAL

Physical model

In this study, the batch protein metal ion binding model in the presence of grafted PEG is based on the isothermal sorption of protein to spherical adsorbent porous particles with an average radius a and a porosity (_i . The adsorption process was conducted in a well-stirred tank with a total system volume V. The liquid volume external to the adsorbent matrix is (_b V and the adsorbent volume is (1-(_b )V. The initial and the transient solute concentration in the liquid are C ₀ and C, respectively. The protein concentrations in the fluid and on the internal solid phase of the adsorbent pores are C _i and Q _i , respectively.

Protein uptake model

The protein batch adsorption is described by a three-resistance model that includes the transport on the film that surrounds the agarose particles, diffusion inside the particles, and reaction (adsorption) with the adsorption sites, in analogy to previous analyses in the literature (Arve and Liapis, 1987; Horstmann and Chase, 1989; Carlsson et al., 1994; Kempe, 2007). The simplifying assumptions of the model are: (1) the bulk solution in perfectly mixed, (2) uniform adsorbent particle size, (3) symmetry in the θ and ϕ directions, (4) negligible adsorption of protein on the outer surface of the adsorbent particle, (5) the adsorption of protein does not change the adsorbent porosity significantly, (6) isothermal conditions, (7) homogeneous diffusion within the adsorbent pores and the boundary layer, (8) homogeneous adsorption and desorption rate coefficients. The size of the adsorbent is considered as an average size, for the sake of model simplicity; also highly monodisperse agarose beads were used (40±8 µm).

Protein molecular diffusivity was assumed constant since the working protein conditions can still be considered as dilute such that the respective activity coefficient can be set to 1.0, making the molecular diffusivity independent of protein concentration (McBride and Rodgers, 2013).

To describe the change of protein concentration in the bulk liquid outside the adsorbent, a mass balance is applied, leading to Eq. 1, which establishes that the protein concentration in the liquid changes due to transport of the protein in the bulk to the surrounding stagnant film that covers the adsorbent particles and characterized by a film mass transfer coefficient.

$$\frac{dC}{dt}=-\frac{3}{a}\frac{\left(1-{\epsilon }_b\right)}{{\epsilon }_b}{k}_L\left(C-{C_i|}_{R=a}\right)$$ (1)

where C is the protein concentration in the bulk liquid (mg/mL), C _i is the protein concentration in the liquid within the adsorbent pores (mg/mL), R is the radial coordinate (m), t is the time coordinate (s), a is the adsorbent radius (m), k _L is the convective mass transfer coefficient (m/s), and ( _b is the volumetric fraction of liquid in the total system volume (-).

The change of protein in the liquid inside the spherical adsorbent particles is given by the mass balance described in Eq. 2. Here the diffusion within the spherical particle, considering only diffusion in the radial direction and under transient conditions is described and characterized by a constant effective diffusivity, D _e .

$${\epsilon }_i\frac{\partial C_i}{\partial t}+(1-{\epsilon }_i)\frac{\partial Q_i}{\partial t}={\epsilon }_i{D}_e\left(\frac{{\partial }^2{C}_i}{\partial R^2}+\frac{2}{R}\frac{\partial C_i}{\partial R}\right)$$ (2)

where Q _i is the protein concentration on the adsorbent surface (mg/mL), D _e is the effective diffusivity of the protein (m²/s), and ( _i is the volumetric fraction of pores in the adsorbent particle (-). The three-resistance model representation is schematically shown in Fig. 1.

Figure 1
Schematic of the three resistances considered for protein adsorption.

Equilibrium of the adsorption process

The second term on the left hand side of Eq. 2 requires a constitutive equation to describe the adsorption process. Several expressions have been proposed for this equation (Limousin et al., 2007). For the case of protein adsorption, a second-order forward reaction and a first-order backward reaction have been used extensively (Chase, 1984; Arnold et al., 1985; Skidmore et al., 1990), and for this work it appeared to be the most appropriate (Eq. 3).

$$\frac{\partial Q_i}{\partial t}=k_1{C}_i\left(q_m-Q_i\right)-k_{-1}{Q}_i$$ (3)

where k ₁ is the adsorption rate constant (mL/mg·s) and k _-1 is the desorption rate constant (1/s). At equilibrium, Eq. 3 >to the familiar Langmuir isotherm (Eq. 4).

$$Q_i^{*}=\frac{q_m{C}_i^{*}}{K_d+C_i^{*}}$$ (4)

where q _m is the maximum protein capacity (mg/mL), K _d is the desorption equilibrium constant, C^*_i is the protein concentration in the liquid at equilibrium (mg/mL), and Q^*_i is the protein concentration in the adsorbent at equilibrium (mg/mL).

Initial and boundary conditions

In order to evaluate the solution of Eqs. 1 - 3, initial and boundary conditions need to be incorporated. These are obtained under the following considerations, at the beginning of the experiment the liquid inside the spherical particles and the surface of the adsorbent are free of protein, while the liquid outside the adsorbent contains the selected initial concentration. Thus, the initial conditions are given in Eqs. 5 - 7.

$$t=\mathrm{0,}C=C_0$$ (5)

$$t=\mathrm{0,}{C}_i=\mathrm{0,}0\le R\le a$$ (6)

$$t=\mathrm{0,}{Q}_i=\mathrm{0,}0\le R\le a$$ (7)

Considering particle symmetry at the center of the spherical particles the following boundary condition is established and given by Eq. 8.

$$R=\mathrm{0,}\frac{\partial C_i}{\partial R}=\mathrm{0,}t>0$$ (8)

A mass balance on the surface of the spherical particles (here, it is assumed that no accumulation occurs or it is negligible) establishes that the mass transfer rate by convection from the bulk fluid to the surface of the particle is equal to the Fickian diffusion mass transfer rate from the surface to the interior of the particle. Thus, the second boundary condition established at the particle surface is given by Eq. 9.

$$R=a,k_L\left(C-C_i\right)=D_e\frac{\partial C_i}{\partial R},t>0$$ (9)

Dimensionless mathematical model

Introducing the dimensionless variables c = C/C ₀, c _i = C _i /C ₀, q _i = Q _i /q _m , r = R/a, and τ = D _e t/a ² to Eqs. 1to3 and Eqs. 5to9 generates the three-resistance model used in this work (Eq. 10 to 20).

$$\frac{\partial c}{\partial \tau }=-A\left(c-{c_i|}_{r=1}\right)$$ (10)

$$\frac{\partial c_i}{\partial \tau }=\frac{{\partial }^2{c}_i}{\partial r^2}+\frac{2}{r}\frac{\partial c_i}{\partial r}-B\frac{\partial q_i}{\partial \tau }$$ (11)

$$\frac{\partial q_i}{\partial \tau }=E{c}_i\left(1-q_i\right)-k_1{K}_d{q}_i$$ (12)

$$\tau =\mathrm{0,}{c}_i=\mathrm{0,}{q}_i=\mathrm{0,}0\le r\le 1$$ (13)

$$\tau =\mathrm{0,}c=1$$ (14)

$$r=\mathrm{0,}{\frac{\partial c_i}{\partial r}|}_{r=0}=\mathrm{0,}\tau >0$$ (15)

$$r=\mathrm{1,}\left(c-{c_i|}_{r=1}\right)=\frac{1}{{\text{Βi}}_{\text{m}}}{\frac{\partial c_i}{\partial r}|}_{r=0},\tau >0$$ (16)

where

$$A=\frac{3\left(1-{\epsilon }_b\right)}{{\epsilon }_b}\frac{k_{L}a}{D_e}$$ (17)

$$B=\frac{q_m}{C_0}\frac{\left(1-{\epsilon }_i\right)}{{\epsilon }_i}$$ (18)

$$E=\frac{k_1{a}^2{C}_0}{D_e}$$ (19)

$${\text{Bi}}_{\text{m}}=\frac{k_{L}a}{{\epsilon }_i{D}_e}$$ (20)

MATERIALS AND METHODS

Adsorbents

The studied adsorbents were based on agarose beads from Inovata. The bead structure is 9.5% cross-linked agarose with a mean particle size of 40 µm. Two types of adsorbents were prepared following protocols described elsewhere (Gonzalez-Ortega, 2010). The adsorption sites were either ion exchange groups or immobilized metal ions. The adsorbents studied are described in Table 1.

Table 1
Adsorbents studied. Super-indices L and H stand for low and high PEG grafting density, respectively.

Adsorbent	Description
IDA-Cu(II)	Grafted iminodiacetic acid (IDA) with Cu(II)
PEG5-IDA-Cu(II)L	Grafted PEG of 5 kDa and IDA with Cu(II)
PEG5-IDA-Cu(II)H	Grafted PEG of 5 kDa and IDA with Cu(II)
Q	Grafted trimethyl ammonium (TMA)
PEG5-QL	Grafted PEG of 5 kDa and TMA
PEG5-QH	Grafted PEG of 5 kDa and TMA
PEG2-QL	Grafted PEG of 2 kDa and TMA
PEG2-QH	Grafted PEG of 2 kDa and TMA

Batch protein adsorption studies

Batch uptake of proteins by the adsorbents studied was performed using 0.05 g of equilibrated gel and 10 mL of a specified protein solution at a given initial concentration. The adsorbents based on immobilized metal ions were previously equilibrated (after saturation with a 20 mM CuSO₄ solution) with a buffer consisting of 20 mM phosphate and 0.5 M NaCl (pH 7.0), while the adsorbents based on ion exchange groups were previously equilibrated (after saturation with a 1 M NaCl solution) with a buffer consisting of 50 mM Tris (pH 8.3).

The adsorption of proteins was performed with agitation at 25°C and liquid samples were withdrawn for absorbance measurements at 280 nm to construct uptake characteristic curves. The experimental data was also used to obtain adsorption isotherms once the equilibrium was reached. To determine the protein concentration adsorbed at the surface of the adsorbent at equilibrium, a mass balance given by Eq. 21 was used.

$$Q_i^{*}=\frac{V_l\rho \left(C_0-C_i^{*}\right)}{m}$$ (21)

where V _l is the liquid volume (mL), ( is the apparent adsorbent density (mg/mL), and m is the wet mass of adsorbent (g).

Model solution and parameter estimation

Eq. 10-12 along with the initial and boundary conditions were simultaneously solved by the method of lines (Davis, 1984). The radial coordinate of the spherical particles was discretized using fourth-order compact finite differences (Lele, 1992). Integration in time of the discretized equations was performed first by the stiff integrator ode15s from MATLAB^® and, once the parameters were close to their optimized values, the integrator was switched to ode45 to increase the global accuracy. Parameter estimation (k ₁ and D _e ) was performed using an unconstrained non-linear optimization with fminsearch from MATLAB^®. The numerical solution scheme for integration and parameter estimation in a MATLAB® environment is shown in Fig. 2. fminsearch is based on a direct search method (Nelder-Mead method) looking to minimize Σ(y _i -y _p )² where y _i and y _p are the experimental and predicted values for c. For the estimation procedure six batch runs were used for each analyzed adsorbent. Since fminsearch does not evaluate confidence intervals for the estimated parameters, each standard error was computed using [Σ(y _i -y _p )²/DF∙cov(i,i)]^0.5 where DF are the degrees of freedom and cov(i,i) is the i-th element of the covariance matrix that was constructed using second-order finite differences. Finally each confidence interval was calculated using σ = 0.05.

Figure 2
Scheme for model solution and parameter estimation.

Internal and external porosities, ( _i and ( _b , were determined using acetone and blue dextran as tracers according to procedures described elsewhere (Lameloise and Viard, 1994). The parameters, q _m and K _d , from the Langmuir isotherms were obtained by non-linear fitting of the experimental equilibrium data. The film mass transfer coefficient, k _L , was calculated using the correlation given in Eq. 22 (Calderbank and Moo-Young, 1961; Geankoplis, 1983).

$$k_L=\frac{D_{AB}}{a}+0.31{\left(\frac{{\mu }_L}{{\rho }_L{D}_{AB}}\right)}^{-2/3}{\left[\frac{\left(\rho -{\rho }_L\right){\mu }_{L}g}{{\rho }_L^2}\right]}^{1/3}$$ (22)

where D _AB is the molecular diffusivity (m²/s), µ _L is the viscosity of the solution (kg/m∙s), and ρ _L is the density of the solution (kg/m³). Eq. 22 is generally accepted to predict the film mass transfer coefficient when using adsorbents with radius smaller than 300 µm.

Implementation of compact finite differences

As previously mentioned, the derivatives in the r-direction were estimated using fourth-order compact finite differences. For this, consider that the adsorbent particle is divided in N equidistant nodes in the r-direction such that the mesh shown in Fig. 3 is generated. For the sake of simplicity N is set to 7 (for the model solution it was set to 100). From Fig. 2, the protein concentration in the liquid within the pores is represented by c with sub-indices 1 to 7, the concentration of protein in the bulk solution is represented as c ₈, and the protein concentration in the surface of the adsorbent is represented by c with sub-indices 9 to 15.

To discretize a first derivative Eq. 23 was used, while Eq. 24 was employed for a second derivative. Both represent fourth-order approximations and h represents the distance between nodes.

$$\frac{1}{4}{c}_{i-1}^{\text{'}}+c_i^{\text{'}}+\frac{1}{4}{c}_{i+1}^{\text{'}}=\frac{3}{4h}\left(c_{i+1}-c_{i-1}\right)$$ (23)

$$\frac{1}{10}{c}_{i-1}^{\text{'}\text{'}}+c_i^{\text{'}\text{'}}+\frac{1}{10}{c}_{i+1}^{\text{'}\text{'}}=\frac{6}{5h^2}\left(c_{i+1}-2c_i+c_{i-1}\right)$$ (24)

To evaluate a first derivative at positions i = 2 to 6 Eq. 25 can be used, while Eq. 26 can be used to evaluate a second derivative. Matrices A1 and A2 and vectors E1 and E2 are determined according to the Appendix and these include the effect of boundary conditions (Eqs. 15 and 16).

$$c_r^{\text{'}}=A1c_r+c_{8}E1$$ (25)

$$c_r^{\text{'}\text{'}}=A2c_r+c_{8}E2$$ (26)

To generate a set of ordinary differential equations (ODEs) at each internal node on the radial mesh; Eq. 11 must be applied. Starting with an example at node i = 2, the next ODE will be generated

$$\frac{d{c}_2}{dt}=c_2^{\text{'}\text{'}}-\frac{2}{h}{c}_2^{\text{'}}-B\left[E{c}_2\left(1-c_{10}\right)-k_1{K}_d{c}_{10}\right]$$ (27)

where C^"₂ can be calculated using Eq. 26 as:

$$c_2^{\text{'}\text{'}}=12.5874c_2-55.1582c_3+61.1407c_4-21.5553c_5+2.9831c_6+0.0022c_8$$ (28)

and C^'₂ using Eq. 25

$$c_2^{\text{'}}=-7.968c_2+11.42c_3-4.2752c_4+0.8759c_5-0.0577c_6+0.005c_8$$ (29)

A similar procedure will generate ODEs for nodes 3 to 6. At node i = 8, from Eq. 10:

$$\frac{d{c}_8}{dt}=-A\left(c_8-c_7\right)$$ (30)

For node i = 9:

$$\frac{d{c}_9}{dt}=E{c}_1\left(1-c_9\right)-k_1{K}_d{c}_9$$ (31)

Equations similar to Eq. 31 will be obtained at nodes 10 to 15.

Figure 3
Example of the mesh for discretization in the r-direction for N = 7.

RESULTS AND DISCUSSION

Adsorption isotherms for adsorbents with immobilized metal ions

The adsorption of proteins was evaluated with the systems IDA-Cu(II) (as control adsorbent), PEG5-IDA-Cu(II)^L (low polymer grafting density), and PEG5-IDA-Cu(II)^H (high polymer grafting density). The model proteins were insulin (INS), ribonuclease A (RNAse A), and bovine serum albumin (BSA), since these are known to have affinity towards immobilized Cu(II) ions with IDA as chelating agent. BSA contains 2-3 histidyl residues on its surface that could interact with immobilized Cu(II) ions, through coordination bonds; also its His 3 residue located at the N-terminus has been recognized as the site for Cu(II) binding (Sulkowski, 1985; Peters and Blumenstock, 1967). RNAse A contains 3 histidyl residues on its surface with Cu(II) chelating properties at positions 12, 105, and 109 (Sulkowski, 1985). INS contains 2 histidyl residues (His 5 and His 10) that could interact with immobilized Cu(II) ions (Farinas et al., 2003).

Fig. 4 shows the isotherms for the studied adsorbents describing the equilibrium adsorption of the three model proteins used in this study. The thermodynamic parameters obtained with experimental and theoretical analyses and the maximum protein capacities and dissociation constants for each case are summarized in Table 2. Fig. 4A describes the adsorption isotherm for BSA with the control chelating adsorbent IDA-Cu(II) where a high affinity between BSA molecules and immobilized Cu(II) ions is evident at low protein equilibrium liquid concentrations. The BSA adsorption isotherms for the adsorbents with grafted PEG were not determined since the adsorption process was forbidden.

Fig. 4B and Fig. 4C describe the adsorption isotherms for RNAse A and INS, respectively. For both proteins favorable behavior for protein adsorption is clearly observed using IDA-Cu(II). The almost linear behavior obtained for INS reveals a high protein capacity for this adsorbent. Fig. 4B also shows the isotherm obtained for RNAse A adsorption on the hybrid system PEG5-IDA-Cu(II)^L; where the presence of a low PEG grafting density reduces the maximum RNAse A capacity by 32% (compared to the system IDA-Cu(II)). RNAse A isotherm with the hybrid system PEG5-IDA-Cu(II)^H was not determined since the high PEG grafting density completely eliminated protein adsorption. Fig. 4C also shows the INS isotherms for the systems PEG5-IDA-Cu(II)^L and PEG5-IDA-Cu(II)^H, where it is evident that protein adsorption is fundamentally affected by the presence of grafted PEG; however, adsorption of INS continues to occur even at a high PEG grafting density.

Table 2
Equilibrium parameters for the proteins studied on adsorbents with immobilized metal ions.

System	BSA		RNAse A		INS
	qm , mg/mL	Kd , mg/mL	qm , mg/mL	Kd , mg/mL	qm , mg/mL	Kd , mg/mL
IDA-Cu(II)	31.699 ±1.809	0.018 ±0.004	26.004 ±9.332	0.390 ±0.295	411.591 ±14.800	0.295 ±0.021
PEG5-IDA-Cu(II)L	nr	nr	17.774 ±3.314	0.244 ±0.120	78.085 ±9.259	0.147 ±0.039
PEG5-IDA-Cu(II)H	nr	nr	nr	nr	18.003 ±2.713	0.011 ±0.008

• nr: not reported since adsorption was negligible

Figure 4
Adsorption isotherms for adsorbents with immobilized metal ions. (A) BSA, (B) RNAse A, (C) INS. (●) IDA-Cu(II), (♦) PEG5-IDA-Cu(II)^L, (▲) PEG5-IDA-Cu(II)^H.

Adsorption isotherms for adsorbents with ion exchange groups

Protein adsorption isotherms were evaluated with the systems Q (as a control adsorbent), PEG5-Q^L (high molecular weight polymer, low polymer grafting density), PEG2-Q^L (low molecular weight polymer, low polymer grafting density), PEG5-Q^H (high molecular weight polymer, high polymer grafting density), and PEG2-Q^H (low molecular weight polymer, high polymer grafting density). The model proteins were BSA, α-lactalbumin (LAC), and INS, since these were expected to interact with the TMA groups on the adsorbents.

Figure 5
Adsorption isotherms for adsorbents with ion exchange groups. (A) BSA, (B) LAC, (C) INS. (●) Q, (♦) PEG5-Q^L, (▲) PEG2-Q^L, (■) PEG5-Q^H.

Fig. 5 shows adsorption isotherms for BSA, LAC and INS for the ion exchange adsorbents. These three proteins are negatively charged under the adsorption conditions and, thus, they are retained by the system Q that contains only anion exchange groups (the quaternary ammonium ions) independent of their molecular weights. At pH 8.3, the amino acids glutamic and aspartic acid are fully deprotonated and, thus, will contribute to the retention of proteins by ion exchange. The proteins studied are retained on the ion exchangers due to the large number of glutamic and aspartic acid residues on their surface (their isoelectric points are below 7.0).

Fig. 5A describes the adsorption isotherm for BSA with the control IEX adsorbent Q where a high electrostatic interaction between the TMA groups and the acid residues from the protein is evident at low protein liquid concentrations at equilibrium. Experiments using the hydrid IEX adsorbents PEG5-Q^L and PEG2-Q^L were performed, but the presence of grafted PEG prevented the interaction of the TMA groups with BSA molecules; thus, an adsorption isotherm could not be determined. A similar situation was obtained when working with systems PEG5-Q^H and PEG2-Q^H. Fig. 5B shows the adsorption isotherms for the adsorbents Q and PEG5-Q^L using LAC. Again when working with the IEX adsorbent Q, LAC molecules strongly interact with the ion exchange groups on the adsorbent surface with the subsequent LAC adsorption even at low protein concentration. When PEG is present in the adsorbent (PEG5-Q^L), LAC continues to adsorb on the hybrid adsorbent, but maximum protein capacity is almost reduced by 90% when compared with the plain IEX adsorbent Q. When working with the hybrid adsorbents PEG5-Q^H, PEG2-Q^L, and PEG2-Q^H; LAC adsorption was almost inexistent and the adsorption isotherms were not determined. Fig. 5C presents adsorption isotherms for INS on the four studied adsorbents. In general the presence of grafted PEG, irrespective of its size or density, reduced the ability of the INS molecules to interact with the TMA groups at equilibrium conditions reducing the maximum protein capacity while increasing the dissociation constant. Despite this, INS continues to adsorb even at a high PEG grafting density only when PEG with a molecular weight of 5 kDa is used. Complete rejection of all studied proteins was obtained with the system PEG2-Q^H. The thermodynamic parameters obtained with experimental and theoretical analyses are summarized in Table 3, where the maximum protein capacities are given for proteins on adsorbents lacking PEG moieties.

Table 3
Equilibrium parameters for BSA, LAC, and INS on adsorbents with ion exchange groups.

System	BSA		LAC		INS
	qm , mg/mL	Kd , mg/mL	qm , mg/mL	Kd , mg/mL	qm , mg/mL	Kd , mg/mL
Q	24.202±2.494	0.098±0.058	31.562±4.616	0.010±0.008	42.332±4.692	0.010±0.007
PEG5-QL	nr	nr	3.628±0.221	0.043±0.012	15.230±1.596	0.052±0.023
PEG5-QH	nr	nr	nr	nr	23.942±8.202	0.693±0.282
PEG2-QL	nr	nr	nr	nr	33.021±13.696	0.808±0.482
PEG2-QH	nr	nr	nr	nr	nr	nr

• nr: not reported since adsorption was negligible

Kinetic studies using adsorbents with immobilized metal ions

BSA kinetics of adsorption

The results for the batch adsorption kinetics of BSA on the adsorbents with immobilized ions are shown in Fig. 6. Fig. 6 compares the behavior of BSA uptake on the plain chelating adsorbent IDA-Cu(II) with the hybrid PEG5-IDA-Cu(II)^L, the low PEG grafting density adsorbent. It is quite evident that the bulk protein concentration drops consistently on the chelating resin, without grafted PEG. Several other initial bulk protein concentrations were tested with similar results (data not shown). These results show very clearly how the protein (BSA) uptake dramatically decreases (no change in the bulk concentration) in the presence of grafted PEG on the hybrid agarose adsorbents.

The continuous line in Fig. 6 represents the solution of the theoretical model (Eq. 10 - 12) after parameter estimation. For the system IDA-Cu(II), an average effective diffusivity (D _e ) of 2.01±0.14(10^-11 m²/s was determined, while the adsorption rate constant (k ₁) provided an average value of 0.72±0.05(10^-3 mL/mg(s. Membrez et al. (1996) experimentally found an effective diffusivity of 2.2(10^-11 m²/s for HSA (human serum albumin) adsorption on similar agarose gels (Sepharose CL-6B). The molecular diffusivity of BSA at 25 ˚C has been reported with a value of 5.94(10^-11 m²/s (Tyn and Gusek, 1990) indicating that the agarose matrix in fact affects protein diffusion considerably. Table 3 shows values of the forward rate constant for BSA and HSA adsorption on several chromatographic systems. Despite the difference in these adsorption systems, the forward rate constant obtained in this work is comparable to values reported in the literature and presented in Table 4.

Figure 6
BSA adsorption kinetics using adsorbents with immobilized metal ions. Open symbols for PEG5-IDA-Cu(II)^L and closed symbols for IDA-Cu(II).

Table 4
Adsorption rate constants for albumin on different chromatographic media.

Adsorbent	Interaction	Protein	k1 , 10-3 mL/mg(s	Reference
Sepharose CL-6B	Dye-Ligand	BSA	1.2	Chase (1984)
Sepharose CL-4B	Dye-Ligand	BSA	3.37-6.12	Arnold and Blanch (1986)
Superose CL-6B	Dye-Ligand	BSA	1.45-1.97	Arnold and Blanch (1986)
Sepharose CL-6B	Dye-Ligand	BSA	0.10-1.58	Boyer and Hsu (1992)
DEAE Cellulose	Anion exchange	BSA	0.3-0.9	Lopez et al. (1997)
Agarose CL 6%	Dye-Ligand	BSA	0.9-2.68	Noriega et al. (1997)
DEAE Sepharose	Anion exchange	HSA	0.49	Aboudzadeh et al. (2006)
IDA-Cu(II)	Chelation	BSA	0.72±0.05	This work
Q	Anion exchange	BSA	0.08±0.007-0.88±0.058	This work

RNAse A kinetics of adsorption

Fig. 7 shows the results and comparison for the batch adsorption of RNAse A on the chelating and hybrid systems. Fig. 7 compares the behavior of RNAse A uptake on the plain chelating adsorbent IDA-Cu(II), the hybrid PEG5-IDA-Cu(II)^L with low PEG grafting density, and the hybrid PEG5-IDA-Cu(II)^H with high PEG grafting density for a selected initial RNAse A concentration. Several other initial concentrations were tested with similar results (data not shown). The presence of grafted PEG clearly affects the adsorption kinetics, reducing the equilibrium protein concentration on the adsorbent. At higher PEG grafting density, as clearly seen in Fig. 7, RNAse A was almost entirely excluded from the immobilized Cu(II) ions and the protein uptake was almost insignificant.

Again, as in the case with BSA adsorption, the continuous lines in Fig. 7 describe the solution of the model (Eq. 10-12) eqparameter estimation. For the plain chelating gel system, Nov-IDA-Cu(II), an average effective diffusivity (D _e ) of 6.93±1.02(10^-11 m²/s was obtained while the adsorption rate constant (k ₁) was determined with an average value of 3.32±0.39(10^-3 mL/mg(s. For the case of the hybrid system PEG5-IDA-Cu(II)^L, average values of 4.91(10^-11±0.14 m²/s and 1.70±0.10(10^-3 mL/mg(s were estimated for the effective diffusivity and the adsorption rate constant; respectively. The molecular diffusivity of RNAse A at 25 ˚C has been reported in the literature with a value of 1.21(10^-10 m²/s (Creeth, 1958) and, as expected, the effective diffusivities obtained in the present work are lower than this reported value. The effect of the presence of grafted PEG on the adsorption of RNAse A is to reduce both the average value of its effective diffusivity and its adsorption rate constant (an almost 50% decrease).

Figure 7
RNAse A adsorption kinetics using adsorbents with immobilized metal ions. (●) PEG5-IDA-Cu(II)^H, (♦) PEG5-IDA-Cu(II)^L, (▲) IDA-Cu(II).

INS kinetics of adsorption

Fig. 8 shows INS batch adsorption kinetics on the chelating and hybrid systems. Fig. 8 describes the INS batch uptake on the chelating system IDA-Cu(II) at a selected initial protein concentration. In this case, as expected and similarly to the other proteins, the bulk protein concentration is greatly reduced, implying that in fact INS has strong affinity for the chelating resin. Fig. 8 also shows the INS adsorption kinetics with the systems having two densities of grafted PEG. Here it is clearly seen that INS continues to adsorb on the immobilized Cu(II) ions although to different extents, however, inversely proportional to the PEG grafting density.

In Fig. 8 the solution of the mathematical model (Eq. 10-12) after parameter estimation, similarly to the previous cases with BSA and RNAse A, is given by the continuous lines. For INS with the system IDA-Cu(II), an average effective diffusivity (D _e ) value of 1.333±0.071(10^-10 m²/s was determined, whereas the adsorption rate constant (k ₁) varied from 0.16 to 0.27±0.02 mL/mg(s, showing a slight dependence on the INS initial liquid concentration. For the case of the hybrid systems with grafted PEG, the effective diffusivity changed only slightly between 1.410±0.064(10^-10 m²/s and 1.328±0.057(10^-10 m²/s for the low and high PEG grafting density, respectively. In the literature, the molecular diffusivity of INS has been reported as 1.63(10^-10 m²/s (Sober et al., 1968) and it appears to be only slightly different; however, no significant changes appear to be associated with the diffusion of the protein with or without grafted PEG. These practically constant values imply that, due to its relatively small size and for the working PEG densities, the permeation effect is negligible in diminishing its diffusion and consequent effective transport to the surface of the adsorbent.

The most significant change observed in the systems for INS adsorption was in the values obtained for the adsorption rate constants for the plain chelating adsorbent and hybrid PEG grafted systems. The adsorption rate constants, in the case of the hybrid systems, were determined with values 2 to 3 orders of magnitude lower than the values obtained with the chelating gels without grafted PEG. These results suggest that, although the PEG polymer does not restrict the transport of the protein to the adsorption sites, the PEG molecules create steric hindrance and thus significantly reduce the effective interaction of the histidyl residues on the protein with the immobilized Cu(II) ions on the adsorbent surface.

Figure 8
INS adsorption kinetics using adsorbents with immobilized metal ions. (●) PEG5-IDA-Cu(II)^H, (♦) PEG5-IDA-Cu(II)^L, (▲) IDA-Cu(II).

Kinetic studies using adsorbents with ion exchange groups

BSA kinetics of adsorption

Fig. 9 shows batch adsorption kinetics for BSA on IEX systems for different initial protein concentrations. In Fig. 9 a comparison is made between BSA uptake by the plain adsorbent Q (solid symbols) and the adsorbent with grafted PEG named PEG5-Q^L (open symbols). For the former a reduction in protein concentration in the liquid is present for all initial protein concentrations, while the latter had no appreciable reduction. Thus, the presence of grafted PEG, even at a low PEG grafting density, completely rejected BSA, forbidding the interaction with the ion exchange adsorption sites. Thus, BSA maximum capacity is reduced from 24.202 mg/mL for the system Q to a level of adsorption that is negligible under the current experimental conditions for the hybrid system PEG5-Q^L. The large size of BSA, compared to the size of grafted PEG, favors its rejection from the polymer volume as described elsewhere (Abbott et al., 1991). The systems PEG2-Q^L, PEG2-Q^H, and PEG5-Q^H were also tested to determine adsorption kinetics, but the results (data not shown) were comparable to the system PEG5-Q^L.

The continuous line in Fig. 5 represents the solution of the three-resistance model (Eq. 10-12) after parameter estimation. The effective diffusivity (D _e ) for BSA was estimated to have an average value of 2.13±0.24(10^-11 m²/s while the adsorption rate constant (k ₁) varies from 0.8±0.07(10^-4 to 8.8±0.58(10^-4 mL/mg(s, increasing as the initial liquid concentration of protein decreases. Table 4 shows values of the forward rate constant for BSA and HSA (human serum albumin) adsorption on several chromatographic systems.

Figure 9
BSA adsorption kinetics using adsorbents with ion exchange groups. Open symbols for PEG5-Q^L and closed symbols for Q.

LAC kinetics of adsorption

The results from batch uptake of LAC using the IEX systems, with and without grafted PEG, are presented in Fig. 10 for a selected initial protein concentration (other LAC initial concentrations were tested with similar results, data not shown). Fig. 10 shows the adsorption kinetics for the plain system Q, where it is clear that the bulk protein concentration is significantly reduced for the selected initial protein concentration. The same Fig. 10 presents the batch uptake results using the hybrid gel PEG5-Q^L (having a low grafting PEG density), where the presence of grafting PEG is clearly reducing the amount of protein that was bound by the hybrid gel. Fig. 10 also shows the results using the high grafting PEG density gel (PEG5-Q^H), where it is clear that the presence of grafted PEG inhibits, almost completely, the adsorption of LAC under the studied experimental conditions. Experiments performed with the systems PEG2-Q^L and PEG2-Q^H gave similar results (data not shown) to the system PEG5-Q^H.

The continuous lines in Figure 10 represent the solution of the model (Eq. 10-12) after parameter estimation. The major dimension of LAC is lower than the size of the PEG (5 kDa) and thus the protein is being somewhat retained by the single PEG modified IEX system, but an increase in PEG density “completely” rejects the protein from the adsorption sites. This somehow indicates that the polymer volume is not being penetrated by the protein and that adsorption in the single modified system could occur on areas of the adsorbent that have a local low PEG grafting density. For the case of the grafted PEG molecule with a molecular weight of 2 kDa, the major dimension of the LAC molecule is now larger and the polymer volume is not being penetrated by the protein with the subsequent result of no appreciable adsorption.

The estimated average effective diffusivities (D _e ) were 0.95±0.07(10^-10 m²/s and 0.84±0.04(10^-10 m²/s for the systems Q and PEG5-Q^L, respectively. The molecular diffusivity for LAC is reported in the literature as 1.06(10^-10 m²/s (Tyn and Gusek, 1990) and it is slightly reduced by the agarose matrix (which is highly porous) and the grafted PEG. The adsorption rate constant (k ₁) was reduced by one order of magnitude (from an average value of 0.0632±0.0102 to 0.0076±0.0018 mL/mg·s) due to the presence of grafted PEG and did not show a dependence on the initial LAC liquid concentration.

Figure 10
LAC adsorption kinetics using adsorbents with ion exchange groups. (●) PEG5-Q^H, (♦) PEG5-Q^L, (▲) Q.

INS kinetics of adsorption

The adsorption kinetics of INS on IEX systems for the same initial protein concentration are shown in Fig. 11. Other initial INS concentrations were also tested with similar results (data not shown). For the plain system Q, it is clear that a significant reduction of bulk protein concentration was obtained due to adsorption by the ion exchange sites. Fig. 11 shows how the presence of grafted PEG reduces the adsorption of INS on hybrid IEX-PEG gels irrespective of grafting density or PEG molecular weight.

The continuous lines in Fig. 11 represent the solution of the model (Eq. 10-12) after parameter estimation. The estimated parameters for all prepared adsorbents are shown in Table 5.

In the literature, the molecular diffusivity of INS seems to be only slightly reduced on the adsorbent systems (with or without grafted PEG), which could be associated with the small size of the INS. The most significant change in the systems with INS adsorption was in the values obtained for the adsorption rate constant for the plain and PEG grafted gels. When compared with the value obtained for the plain adsorbent Q, the systems having grafted PEG presented values with a decrease of one order of magnitude. This result suggests that, although the PEG does not restrict the transport of the protein to the adsorption sites, the PEG molecule significantly affects the effective interaction between the protein and the TMA groups.

Figure 11
INS adsorption kinetics using adsorbents with ion exchange groups. (●) PEG5-Q^H, ((■) PEG2-Q^L, (♦) PEG5-Q^L, (▲) Q.

Table 5
Optimized parameters for INS adsorption kinetics on adsorbents with ion exchange groups.

System	k1 , 10-3 mL/mg(s	De , 10-10 m2 /s
Q	22.93±2.37	1.31±0.09
PEG5-QL	4.72±1.24	1.26±0.04
PEG2-QL	1.23±0.12	1.04±0.17
PEG5-QH	1.35±0.28	0.97±0.08

Error analyses

Two error analyses were performed to determine the order of convergence (OC). Since both the discretization in the r-direction and the integration in time used fourth-order approximation schemes, the error analyses should predict an order of convergence around four. The first error analysis was conducted calculating the relative error (RE) between the numerical solution (obtained with different values of h) and a numerical solution for N = 500 (considered as an analytical solution) for c _i at r = 1 using different grid sizes. r = 1 was selected since larger concentration gradients are present at the pore entrance. Table 6 shows the results of this analysis; the order of convergence was calculated with Eq. 32.

$$OC=\frac{\log \left({\text{RE}}_{\text{2}}/{\text{RE}}_{\text{1}}\right)}{\log \left(h_2/h_1\right)}$$ (32)

where RE₁ and RE₂ are the absolute errors corresponding to grid sizes h ₁ and h ₂. For the second error analysis the protein adsorption kinetics were considered. Eq. 32 was still used to calculate the order of convergence using the summation of the relative errors. Table 6 also shows the results of the second error analysis corroborating the use of fourth-order schemes.

Table 6
Error analyses to determine the order of convergence.

N	h	First analysis		Second analysis
		RE	OC	RE	OC
7	0.167	2.59×10-5	-	7.90×10-4	-
9	0.125	8.28×10-6	3.96	2.57×10-4	3.91
11	0.100	3.12×10-6	4.37	1.10×10-4	3.79
13	0.083	1.55×10-6	3.85	5.72×10-5	3.61

CONCLUSIONS

A three-resistance model was used to successfully describe the batch adsorption of proteins using regular adsorbents (plain IDA-Cu(II) and Q) and hybrid adsorbents bearing grafted PEG recognized as restricted access media. The effective diffusivity and the adsorption rate constant as parameters of the model were sufficient to describe the adsorption kinetics of several proteins with molecular weights in the range 5.8-67.0 kDa. The model allowed discriminating the controlling process behind the restricted access condition created by the grafted PEG molecules, which is associated with the interaction between the protein molecule and the adsorption site that reduces or inhibits the adsorption process. Thus, the model predicts a reduction in the adsorption rate constant that is even noticeable for small molecular weight grafted PEG or a high PEG grafting density. The use of compact finite differences allowed obtaining a highly accurate solution using a three-point stencil that accelerates the solution and the total time for the computation of the parameters being optimized.

NOMENCLATURE

• a  - Radius of adsorbent particle, m
• c  - Dimensionless protein concentration in the bulk liquid
• C  - Protein concentration in the bulk liquid, mg/mL
• c _i  - Dimensionless protein concentration in the liquid within the pores
• C _i  - Protein concentration in the liquid within the pores, mg/mL
• C _i ^*  - Protein concentration in the liquid within the pores at equilibrium, mg/mL
• C ₀  - Initial protein concentration in the bulk liquid, mg/mL
• D _e  - Effective protein diffusivity, m²/s
• k _L  - Convective mass transfer coefficient, m/s
• k ₁  - Adsorption rate constant, mL/mg·s
• k  _-1 - Desorption rate constant, 1/s
• K _d  - Desorption equilibrium constant, mg/mL
• q _i  - Dimensionless protein concentration on the adsorbent
• q _m  - Maximum protein capacity, mg/mL
• Q _i  - Protein concentration on the adsorbent, mg/mL
• Q _i ^*  - Protein concentration on the adsorbent at equilibrium, mg/mL
• r  - Dimensionless radial coordinate -
• R  - Radial coordinate, m
• Bi_m  - Mass transfer Biot number
• t  - Time coordinate, s
• V  - Total system volume, m³
Greek letters
• ε _b  - Volumetric fraction of liquid in the total system volume
• ε _i  - Volumetric fraction of pores in the adsorbent particle
• (  - Dimensionless time coordinate

REFERENCES

• Abbott, N.L., Blankschtein, D. and Hatton, T.A. Protein Partitioning in Two-Phase Aqueous Polymer Systems. 1. Novel Physical Pictures and a Scaling-Thermodynamic Formulation. Macromolecules, 24, 4334 (1991).
• Arnold, F.H., Blanch, H.W. and Wilke, C.R., Analysis of Affinity Separations I: Predicting the Performance of Affinity Adsorbers, The Chemical Engineering Journal, 30, B9-B23 (1985).
• Arve, B.H. and Liapis, A.I., Modeling and Analysis of Biospecific Adsorption in a Finite Bath, AIChE Journal, 33, 179 (1987).
• Calderbank, P.H. and Moo-Young, M.B., The Continuous Phase Heat and Mass-Transfer Properties of Dispersions, Chemical Engineering science, 16, 39 (1961).
• Carlsson, F., Axelsson, A. and Zacchi, G., Mathematical Modelling and Parametric Studies of Affinity Chromatography, Computers & Chemical Engineering, 18, S657 (1994).
• Chase, H.A., Prediction of the Performance of Preparative Affinity Chromatography, Journal of Chromatography, 297, 179 (1984).
• Creeth, J.M., Studies of Free Diffusion in Liquids with the Rayleigh Method. III. The Analysis of Known Mixtures and Some Preliminary Investigations with Proteins, Journal of Physical Chemistry, 62, 66 (1958).
• Davis, M.E., Numerical Methods and Modeling for Chemical Engineers. John Wiley & Sons, New York (1984).
• Farinas, C.S., Bueno, S.M.A. and Miranda, E.A., Adsorption of Glucagon and Insulin on an Immobilized Metal Ion Affinity Chromatography Silica Matrix, Adsorption Science & Technology, 21, 883 (2003).
• Geankoplis, C.J., Transport Processes and Unit Operations. Allyn & Bacon, New York, (1983).
• Gonzalez-Ortega, O., Study and Characterization of Dual-Function Affinity Chromatographic Adsorbents Having Size Exclusion and Adsorption Properties to Isolate, Purify and Recover Small Biomolecules from Complex Biological Mixtures. Ph.D. Thesis, The University of Arizona (2010).
• González-Ortega, O., Porath, J. and Guzmán, R., Adsorption of Peptides and Small Proteins with Control Access Polymer Permeation to Affinity Binding Sites. Part I: Polymer Permeation-Immobilized Metal Ion Affinity Chromatography Separation Adsorbents with Polyethylene Glycol and Immobilized Metal Ions, Journal of Chromatography A 1227, 115 (2012a).
• González-Ortega, O., Porath, J. and Guzmán, R., Adsorption of Peptides and Small Proteins with Control Access Polymer Permeation to Affinity Binding Sites. Part II: Polymer Permeation-Ion Exchange Separation Adsorbents with Polyethylene Glycol and Strong Anion Exchange Groups, Journal of Chromatography A, 1227, 126 (2012b).
• Gurarslan, G., Numerical Modelling of Linear and Nonlinear Diffusion Equations by Compact Finite Difference Method, Applied Mathematics and Computation, 216, 2472 (2010).
• Gutierrez, R., Del Valle, E.M.M. and Galan, M.A., Characterization of Mass Transport Processes in IMAC Chromatography by Dynamics Methods, Biochemical Engineering Journal, 35, 264 (2007).
• Horstmann, B.J. and Chase, H.A., Modeling the Affinity Adsorption of Immunoglobulin G to Protein A Immobilised to Agarose Matrices, Chemical Engineering Research & Design, 67, 243 (1989).
• Kempe, H., Modeling and Experimental Determination of Chromatographic Parameters by Frontal Chromatography. Ph.D. Thesis, Stockholm University (2007).
• Lameloise, M.L. and Viard, V., Choice of a Tracer for External Porosity Measurement in Ion-Exchange Resin Beds, Journal of Chromatography A, 679, 255 (1994).
• Lele, S.K., Compact Finite Difference Schemes with Spectral-like Resolution, Journal of Computational Physics, 103, 16 (1992).
• Membrez, J., Infelta, P.P. and Renken, A., Use of the Laplace Transform Technique for Simple Kinetic Parameters Evaluation. Application to the Adsorption of a Protein on Porous Beads, Chemical Engineering Science, 51, 4489 (1996).
• Limousin, G., Gaudet, J.P., Charlet, L., Szenknect, S., Barthes, V. and Krimissa, M., Sorption Isotherms: A Review on Physical Bases, Modeling and Measurement, Applied Geochemistry, 22, 249 (2007).
• McBride, D.W. and Rodgers, V.G., Predicting the Activity Coefficients of Free-Solvent for Concentrated Globular Protein Solutions Using Independently Determined Physical Parameters, PloS one, 8, e81933 (2013).
• Peters, T. and Blumenstock, F.A., Copper-Binding Properties of Bovine Serum Albumin and its Amino-Terminal Peptide Fragment, The Journal of Biological Chemistry, 242, 1574 (1967).
• Sharma, S. and Agarwal, G.P., Adsorption Equilibrium and Kinetics of Egg-White Proteins on Immobilized Metal Ion Affinity Gels for Designing Fractionation, Adsorption, 8, 203 (2002).
• Skidmore, G.L., Horstmann, B.J. and Chase, H.A., Modeling Single-Component Protein Adsorption to the Cation Exchanger S Sepharose FF, Journal of Chromatography , 498, 113 (1990).
• Sober, H.A., Harte, R.A. and Sober, E.K., Handbook of Biochemistry: Selected Data for Molecular Biology. CRC Press, Ohio (1968).
• Sulkowski, E. Purification of proteins by IMAC. Trends in Biotechnology. 3: 1-7 (1985).
• Tyn, M.T. and Gusek, T.W., Prediction of Diffusion Coefficients of Proteins, Biotechnology and Bioengineering, 35, 327 (1990).

Appendix

The procedure to obtain matrices and vectors needed to evaluate first derivatives in the r-direction using a compact finite difference scheme is described. A similar procedure will generate matrices and vectors needed to discretize a second derivative in the r-direction. For discretization of boundary conditions at nodes 1 and 7, for the sake of simplicity, a forward fourth-order finite difference and a backward fourth order finite difference were used to generate Eq. S.1 and S.2, respectively:

$$c_1^{\text{'}}=\frac{-25c_1+48c_2-36c_3+16c_4-3c_5}{12h}=0$$ (S.1)

$$c_7^{\text{'}}=\frac{3c_3-16c_4+36c_5-48c_6+25c_7}{12h}={\text{Bi}}_{\text{m}}\left(c_8-c_7\right)$$ (S.2)

Solving for c ₁ and c ₇:

$$c_1=\frac{48c_2-36c_3+16c_4-3c_5}{25}$$ (S.3)

$$c_7=\frac{12h{\text{Bi}}_{\text{m}}{c}_8-3c_3+16c_4-36c_5+48c_6}{25+12h{\text{Bi}}_{\text{m}}}$$ (S.4)

For node i = 2, applying Eq. 22:

$$\frac{1}{4}{c}_1^{\text{'}}+c_2^{\text{'}}+\frac{1}{4}{c}_3^{\text{'}}=\frac{3}{4h}\left(c_3-c_1\right)$$ (S.5)

Substituting Eq. S.1 and S.3 into Eq. S.5:

$$c_2^{\text{'}}+\frac{1}{4}{c}_3^{\text{'}}=-\frac{144}{100h}{c}_2+\frac{183}{100h}{c}_3-\frac{48}{100h}{c}_4+\frac{9}{100h}{c}_5$$ (S.6)

For node i = 3, applying Eq. 22:

$$\frac{1}{4}{c}_2^{\text{'}}+c_3^{\text{'}}+\frac{1}{4}{c}_4^{\text{'}}=\frac{3}{4h}\left(c_4-c_2\right)$$ (S.7)

Similar equations will be obtained for nodes 4 and 5. For node i = 6, applying Eq. 22:

$$\frac{1}{4}{c}_5^{\text{'}}+c_6^{\text{'}}+\frac{1}{4}{c}_7^{\text{'}}=\frac{3}{4h}\left(c_7-c_5\right)$$ (S.8)

Substituting Eq. S.2 and S.4

$$c_5^{\text{'}}+4c_6^{\text{'}}=-\frac{3F}{h}{c}_3+\frac{16F}{h}{c}_4-\frac{\left(3+36F\right)}{h}{c}_5+\frac{48F}{h}{c}_6+{\text{Bi}}_{\text{m}}\left(12F-1\right)c_8$$ (S.9)

where

$$F=\frac{3{\text{Bi}}_{\text{m}}{}^{-1}+h}{25{\text{Bi}}_{\text{m}}{}^{-1}+12h}$$ (S.10)

In matrix notation we can write

$$A{c}_r^{\text{'}}=B{c}_r+c_{8}e$$ (S.11)

where

$$A=\left[\begin{array}{ccccc}1& 1/4& 0& 0& 0\\ 1/4& 1& 1/4& 0& 0\\ 0& 1/4& 1& 1/4& 0\\ 0& 0& 1/4& 1& 1/4\\ 0& 0& 0& 1& 4\end{array}\right],e=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ {\text{Bi}}_{\text{m}}\left(12F-1\right)\end{array}\right],c_r^{\text{'}}=\left[\begin{array}{c}{c}_2^{\text{'}}\\ c_3^{\text{'}}\\ c_4^{\text{'}}\\ c_5^{\text{'}}\\ c_6^{\text{'}}\end{array}\right],B=\frac{1}{h}\left[\begin{array}{ccccc}-144/100& 183/100& -48/100& 9/100& 0\\ -3/4& 0& 3/4& 0& 0\\ 0& -3/4& 0& 3/4& 0\\ 0& 0& -3/4& 0& 3/4\\ 0& -3F& 16F& -\left(3+36F\right)& 48F\end{array}\right]$$

Multiplying Eq. S.11 by the inverse matrix of A allows solving for vector

$$c_r^{\text{'}}=A^{-1}B{c}_r+c_8{A}^{-1}e=A1c_r+c_{8}E1$$ (S.12)

where (for Bi_m = 31.25 as example):

$$A1=\left[\begin{array}{ccccc}-7.968& 11.420& -4.275& 0.876& -0.058\\ -2.688& -1.760& 5.581& -1.344& 0.231\\ 0.720& -4.380& -0.048& 4.498& -0.865\\ -0.192& 1.281& -5.388& 1.351& 3.230\\ 0.048& -0.742& 3.598& -9.904& 5.947\end{array}\right],E1=\left[\begin{array}{c}0.005\\ -.020\\ 0.076\\ -0.282\\ 1.053\end{array}\right]$$

Publication Dates

• Publication in this collection
  Jan 2018

History

• Received
  30 June 2016
• Reviewed
  21 Oct 2016
• Accepted
  23 Nov 2016