Article pubs.acs.org/IECR
Adsorption of Proteins onto Ion-Exchange Chromatographic Media: A Molecular Dynamics Study Juan Liang,*,† Georg Fieg,† Frerich J. Keil,‡ and Sven Jakobtorweihen‡ †
Institute of Process and Plant Engineering and ‡Institute of Chemical Reaction Engineering, Hamburg University of Technology, 21073 Hamburg, Germany ABSTRACT: Molecular dynamics (MD) simulations with coarse-grained models are used to investigate the adsorption process of proteins onto an ion-exchange chromatographic medium. The adsorption of human serum albumin (HSA) and bovine hemoglobin (bHb) on the anion exchanger Q Sepharose FF is studied. These two proteins have similar molecule sizes with different isoelectric points. To obtain a reliable set of data, simulations with different initial conditions are carried out for each protein. Convincingly, the results of the MD simulations are qualitatively consistent with those of previous experiments. HSA reaches a stable adsorption in all simulations, while the binding sites can differ. In contrast, bHb reaches a stable adsorption only in one simulation. Due to the stronger protein−ligand interaction of HSA, it adsorbs faster, stronger, and with more binding sites onto Q Sepharose FF than bHb. The spatial movements during the adsorption process and the adsorbed states are investigated in detail. Moreover, we studied the underlying physical relevance of two parameters (characteristic charge and steric factor) of the steric mass action (SMA) model for HSA and for bHb with MD simulations and compared them with the results of experiments. The results reveal that MD simulations are useful to interpret the physical consistence of the SMA parameters.
1. INTRODUCTION Ion-exchange chromatography (IEC) is widely used for separation and purification of proteins.1 The target protein adsorbs onto an adsorbent via an ion-exchange mechanism. Separation and purification are realized due to different electrostatic interactions between the solutes and the ionexchange chromatographic media. For an improved protein separation by IEC, it is very important to understand the adsorption of proteins onto ion-exchange chromatographic media in detail. Recently, much attention has been paid to adsorption of proteins, and models were proposed to describe it. The steric mass action (SMA) model is one of the most widely used approaches for ion-exchange adsorption of proteins. It was proposed by Brooks and Cramer,2 taking into account the multipointed nature of protein binding and steric hindrance of salt counterions upon protein binding. For a certain ionexchange adsorbent, whose total capacity is Λ, the SMA model considers ion-exchange adsorption of proteins as a stoichiometric exchange. According to the SMA model, the amount of the binding sites is equal to the characteristic charge (υ). Upon protein binding a monovalent salt counterion is displaced for each binding site; further, it results in the steric hindrance of some other salt counterions. The number of the sterically hindered counterions is equal to the steric factor (σ) of the protein. The SMA equation is given as C=
⎞υ Cs ⎛ Q ⎞⎛ ⎜ ⎟⎜ ⎟ ⎝ K ⎠⎝ Λ − ( υ + σ )Q ⎠
The SMA model has been widely reported for predicting and describing nonlinear ion-exchange adsorption of proteins;3−5 it has been applied to ion-exchange displacement systems6−9 and metal-affinity systems.10,11 The model parameters and accuracy of the SMA model in describing ion-exchange adsorption are influenced by pH and ionic strength,12−15 as well as by the ligand density of the ion-exchange media.16 To obtain the SMA parameters, a simple experimental protocol was presented by Seidel-Morgenstern:17 a set of linear elution experiments and nonlinear frontal experiments for both proteins and displacer were carried out. Barz et al.18 obtained the SMA parameters from static batch experiments and showed that all the SMA parameters can be identified by a relatively small number of experiments. Due to the lack of microscopic experimental techniques, the SMA parameters were determined indirectly from experiments, although the parameters have a physical basis. Molecular dynamics (MD) simulations are a powerful tool to study dynamics on an atomistic scale and, thus, can provide direct microscopic information of ion-exchange chromatographic adsorption of proteins within adsorbent pores. The influence of the mobile phase on liquid chromatographic interactions has been investigated with MD simulations.19−21 Furthermore, MD simulations have been used to understand the interactions between proteins and immobilized ligands,21−24 as well as protein conformation transition in chromatography.21,25,26 All-atom models are used in molecular simulations due to their high precision. However, the molecular simulation
(1)
where K is the adsorption equilibrium constant; Q and C are the protein concentrations on the adsorbent (adsorption density) and in the liquid phase, respectively. Cs is the concentration of salt in the liquid phase. © 2012 American Chemical Society
Received: Revised: Accepted: Published: 16049
May 29, 2012 November 5, 2012 November 16, 2012 November 16, 2012 dx.doi.org/10.1021/ie301407b | Ind. Eng. Chem. Res. 2012, 51, 16049−16058
Industrial & Engineering Chemistry Research
Article
Figure 1. Ligand and ion-exchange system model: (a) sketch of atomistic to coarse-grained mapping of a ligand; (b) snapshot of the ion-exchange system at the beginning of a simulation (protein−ligand minimum distance is set to 1.0 nm). For a better view water particles are not shown here; positively and negatively charged particles are colored with red and blue, respectively. The pink arrows represent six initial orientations of the protein.
methods employing all-atom force fields must be appropriately set up and use parameters suitable for protein adsorption. Latour27 has discussed the advantages and drawbacks of MD simulations for protein adsorption in detail. The main drawbacks of all-atom simulations are the limited time and length scales. A promising step toward larger systems and larger time scales is simulations based on coarse-grained (CG) models. Within these models the number of interactions are reduced by using effective interactions for groups of atoms and, thus, simplifying the microscopic system,28,29 by the cost of losing some atomistic details. However, with CG models it is possible to investigate larger systems on a larger time scale. With these models it is still possible to study the driving forces of ion-exchange chromatographic processes21 (i.e., electrostatic interactions). Though adsorption of proteins onto ion-exchange chromatographic media has been investigated by many researchers, little is known about the dynamics within the adsorbent pores due to the lack of powerful microscopic experimental techniques. Therefore, in this work MD simulations with a CG model (MARTINI force field30) was applied to investigate ionexchange adsorption of proteins onto chromatographic media. The adsorption onto Q Sepharose FF of two different proteins, human serum albumin (HSA) and bovine hemoglobin (bHb), was investigated. Q Sepharose FF is a strong anion exchanger, whose ion-exchange ligand is a quaternary amine group. The proteins investigated are typical representatives for proteins which have similar molecular sizes but different isoelectric points. HSA is composed of one amino acid chain, which is composed of 582 amino acids. bHb is composed of four amino acid chains (chains A, B, C, and D). Chains A and C are hemoglobin alpha chains, and each is composed of 141 amino acids. Chains B and D are hemoglobin beta chains, and 146 amino acids are found in each of these chains. The aim of this work is to investigate the adsorption process of proteins onto ion-exchange chromatographic media on a microscopic scale. To check the influence of initial conditions and obtain a reliable set of data, simulations with different initial conditions are carried out for each protein. The spatial movements during the adsorption process and the adsorbed
states are investigated in detail. Two SMA parameters, the characteristic charge, υ, and the steric factor, σ, of both proteins are obtained. The results of the MD simulations are compared with those of previous experiments.12,31,32
2. METHODS 2.1. Molecular Model. Considering the large size of the ion-exchange systems, the MARTINI coarse-grained force field30,33,34 was used in this work. This force field was introduced as a general coarse-grained force field for biomolecular systems. So far it has been used to study a wide range of systems. In general, a rule of mapping four heavy atoms onto one CG is followed, and hydrogen atoms are neglected because of their small size and mass. The model considers four main types of CG particles: polar (P), nonpolar (N), apolar (C), and charged (Q). The ligand of Q Sepharose FF is a quaternary amine group, which is a strong anion exchanger. As shown in Figure 1a, the 14 heavy atoms of the ligand were modeled with three MARTINI particles, whereas one particle is positively charged. This CG representation of the ligand has been tested by comparing results from CG and atomistic simulations (results not shown). CG models of the proteins are generated according to the protein mapping principle of the MARTINI force field,34 starting from the atomistic structures of the Protein Data Bank (HSA, PDB ID 1AO6, in which chain B, which is identical to chain A, was deleted; bHb, PDB ID 1G0A). Only the protein structures of the PDB files were used; all other atoms in the crystal structure were removed. On condition of this work (pH 7.0, Cs = 230 mM), HSA has a charge of −15 and bHb has a charge of +2. Standard CG water in the MARTINI force field is modeled as P4 particles. In order to prevent the CG water from freezing, 10% of the water particles are modeled as BP4.30 The ion particles follow the standard force field. The simulated ion-exchange system contains a protein, the ion-exchange ligands, and the chromatographic base, as shown in Figure 1b. An infinite flat plate model is used to describe the chromatographic base with ligands immobilized onto it, ignoring the radian of the wall. At the bottom of the simulation box (19.604 nm × 20.404 nm × 20 nm), the chromatographic 16050
dx.doi.org/10.1021/ie301407b | Ind. Eng. Chem. Res. 2012, 51, 16049−16058
Industrial & Engineering Chemistry Research
Article
Figure 2. Potential energy between the protein and the ligands plotted over simulation time for different starting orientations (O1−O6): (a) O1 (black), O2 (blue), and O3 (pink) of HSA; (b) O4 (red), O5 (green), and O6 (gray) of HSA; (c) O1 (black), O2 (blue), and O3 (pink) of bHb; (d) O4 (red), O5 (green), and O6 (gray) of bHb.
protein were simulated.23 The six orientations will be represented in this work as O1−O6 (see Figure 1b). 2.2. Simulation Method. The simulations were carried out with GROMACS 4.5.36 Periodic boundary conditions in the x and y directions were used. In order to use pressure coupling, another wall is set at the other side of the box (opposite to the ligand surface), defined as P4 particle type. The basic parameters for the simulations are set as in previous simulations employing the MARTINI force field.30,34 A leapfrog algorithm for integrating Newton’s equations of motion is used with a simulation time step of 25 fs. This time step size is usually used with the MARTINI model, as it was used during parametrization.33,34 The cutoff distance of the nonbonded interactions are rcut = 1.2 nm. A Lennard-Jones shifting from 0.9 to 1.2 nm was used. As the electrostatic interaction is the main driving force for the ion-exchange system, and as a reaction field37 can be used within the MARTINI force field, this method is used here with a long-range dielectric constant εrf = 7838,39 and a short-range dielectric constant εr = 15.30,34 The reference temperature was coupled (coupling constant 2.0 ps) to T = 298.15 K using a Berendsen thermostat. The pressure was coupled only in the z direction (coupling constant 5.0 ps, compressibility 1 × 10−5 bar−1). In this work, all systems have been energy minimized before starting a MD simulation. MD simulations of ion-exchange adsorption of proteins contained two parts: (1) 100 ns equilibration of chromatographic media, where the positions of the protein particles were restrained; (2) 400 ns MD simulation with an unrestrained protein. Only the results of the second part are discussed. Most of the results are calculated and analyzed by the GROMACS suite of programs. To characterize the potential energy between the protein and the ligands, the Lennard-Jones and Coulomb terms between these two groups are used.
matrix with immobilized ligands can be seen in Figure 1b. A wall is set at z = 0 nm to mimic the matrix of the chromatographic media, and 1116 ligands are put as a uniform layer on top of the wall, with the roots restrained to mimic the immobilization. According to the polarity of Sepharose, the atom type of the wall is defined as P4, following the mapping principle of the MARTINI force field.30 The surface density of the ligands is 2.67/nm2 (which is equal to 1.88 mmol of Cl−/ mL of media, according to the total ionic capacity and the specific surface area of Q Sepharose FF).35 The protein is placed into the liquid phase with a minimum distance to the ligands of 1.0 nm (see Figure 1b). The box is further filled with water particles to a density of 1 g/cm3. A certain number of counterions (1101 chlorides for HSA and 1118 chlorides for bHb) are added to keep the system electrically neutral. The CG simulation of HSA contains a total of 67 513 particles (1297 representing HSA, 3348 ligands, 62 868 water and ions), while in the same box of an all-atom simulation there would be a total of 776 609 atoms (5843 HSA atoms, 18 972 ligand atoms, 751 794 water and ion atoms). That is, by using the MARTINI CG model the number of interaction centers is reduced by more than 90%. The CG system of bHb contains 66 347 particles. For a better assessment of structural changes during adsorption, MD simulations of HSA and bHb solvated by water were also performed with the CG model. Note that the dynamics in CG simulations are faster than all-atom dynamics as the friction is lowered by CG models.34 As time-dependent properties, such as diffusion coefficients, are not studied in this work, we show the simulation time in the figures, rather than trying to map them to real times. The speed-up factor for the system under study is anyway uncertain. Since the protein is initially put near the chromatographic base, the diffusion of the protein toward the wall is not studied. To study the influence of the initial orientation of the protein, six different initial orientations (with 90° rotation) for each 16051
dx.doi.org/10.1021/ie301407b | Ind. Eng. Chem. Res. 2012, 51, 16049−16058
Industrial & Engineering Chemistry Research
Article
Figure 3. Distance of protein center of mass to ligands in the z direction (black), minimum distance between protein and ligands (red), and rotation of protein (secondary y-axis) around x-axis (blue) and y-axis (pink) of HSA adsorption, all plotted over simulation time, with different initial orientations (notation of orientations see Figure 1): (a) O1, (b) O2, (c) O3, (d) O4, (e) O5, and (f) O6.
3. RESULTS AND DISCUSSION 3.1. Protein Adsorption. During the simulations the protein moves to the ligands and adsorbs due to the attractive protein−ligand interaction. Meanwhile, adjustment of the protein orientation is needed both before and after adsorption. Before the adsorption, rotation and adjustment lead to a favorable orientation for adsorption. In order to reach a more stable adsorption, the protein further reorients with respect to the ligand surface after the protein is adsorbed. For a better understanding of the adsorption and the reorientation of the protein, the potential energy between the protein and the ligands of both HSA and bHb were calculated (see Figure 2). The minimum distance between the protein and the ligands, the distance of the protein center of mass and the ligands in the z direction, and the rotation of the protein are shown in Figures 3 (HSA) and 4 (bHb). Since the rotation of the protein around the z-axis does not affect the protein−ligand interaction, only protein rotations around the x- and y-axes are shown and discussed here. Moreover, the individual rotations of every chain of bHb were also analyzed (results not shown); there is no obvious difference compared to the protein rotations as a whole. HSA moves toward the ligands for all different initial orientations (see Figure 3); both the distance of the protein
center of mass in the z direction and the minimum distance between the protein and the ligands is decreasing. For O4 stronger rotation is found compared to the other five initial orientations. Although HSA reaches a stable adsorption (see Figure 2b), the protein−ligand interaction is weaker compared to the other orientations (see Figure 2). Therefore, larger rotations are needed to reach a stronger adsorption. HSA stays close to the ligands in the O6 simulation and reaches adsorption quickly within 10 ns. For O2 a more complex behavior is observed: first the protein moves away from the ligands and then starts to move toward the ligands until 64 ns. This is because, at the beginning, the initial orientation is not favorable for adsorption. The protein moves away and rotates and hence reaches a more favorable configuration. Then, the protein stays near the ligands and reaches a stable adsorption around 176 ns. This can also be seen from Figure 2a: the potential energy of the protein−ligand interaction is already negative at 64 ns, but from 176 ns on, the potential energy decreases indicating a stronger adsorption. For all six different initial orientations, the protein center of mass keeps further fluctuating and gets even closer to the ligands after HSA has been already adsorbed (see Figure 3). As the HSA protein has many negatively charged residues distributed all over its surface, further orientational adjustment (rotation and center of mass 16052
dx.doi.org/10.1021/ie301407b | Ind. Eng. Chem. Res. 2012, 51, 16049−16058
Industrial & Engineering Chemistry Research
Article
Figure 4. Distance of protein center of mass to ligands in the z direction (black), minimum distance between protein and ligands (red), and rotation of protein (secondary y-axis) around the x-axis (blue) and y-axis (pink) of bHb adsorption, all plotted over simulation time, with different initial orientations (notation of orientations see Figure 1): (a) O1, (b) O2, (c) O3, (d) O4, (e) O5, and (f) O6.
between bHb and the ligands keeps fluctuating; the bHb center of mass is getting slightly closer to the ligands for the different initial orientations. For O1, O3, O4, and O5, though the protein−ligand interaction energy is negative, bHb keeps fluctuating in the neighborhood of the ligands without adsorption within 400 ns, indicated by the fluctuating protein−ligand minimum distance and the higher potential energy of the protein−ligand interaction. The interaction energies for these simulations are higher than those for HSA. For O6, after fluctuating in the neighborhood of the ligands without stable adsorption, the protein moves away from the ligands at 210 ns. This shows obviously the weak adsorption of bHb. O2 is the only initial orientation from which bHb reaches adsorption. At the beginning of the simulation the minimum bHb−ligand distance is immediately lower, but bHb has not reached a stable adsorption, indicated by the high interaction energy and by the fluctuating minimum distance. Therefore, initially bHb keeps rotating around both the x- and y-axes in one direction to get a better orientation for adsorption, the bHb center of mass moves slowly toward the ligands, and the interaction energy decreases strongly and stays at a low value. bHb adsorbs after 84 ns. Thereafter, only small rotations are observed for further adjustment. Note that the protein movement is limited once it is adsorbed. The center of mass
movement) can be observed. A stable adsorption is always observed for HSA regardless of the initial protein orientations considered, indicated by a lowered protein−ligand interaction energy (see Figure 2 a,b). After the adsorption of HSA the protein−ligand interaction energy is further decreasing, indicating a more stable adsorption after the adjustment of protein position and orientation. For all six different initial orientations, the rotations around both the x- and y-axes are less than 110° (see Figure 3). Note that HSA (isoelectric point = 4.9)12 is highly negatively charged at pH 7.0 and can strongly adsorb onto anion-exchange media due to strong electrostatic interactions, already observed in experiments.12,14,31,32 The negative charges are distributed all over the HSA surface. Hence, in all orientations there are negatively charged AA residues oriented toward the ligands. Therefore, not much rotation is needed before its adsorption. Once it is adsorbed, HSA is bound strongly onto the ligands with a few binding sites. In consequence, large rotations are hindered. For bHb (see Figures 2 and 4), the simulation results show that adsorption of bHb is influenced by the initial protein orientation. Different protein orientations of bHb lead to different results of its adsorption onto the ion-exchange media. A stable bHb adsorption is only observed in the simulation starting from O2. Except for O6 the minimum distance 16053
dx.doi.org/10.1021/ie301407b | Ind. Eng. Chem. Res. 2012, 51, 16049−16058
Industrial & Engineering Chemistry Research
Article
Figure 5. Minimum distance between every amino acid and the ligands over simulation time from MD simulations of HSA adsorption with different initial orientations (notation of orientations see Figure 1): (a) O1, (b) O2, (c) O3, (d) O4, (e) O5, and (f) O6.
of bHb keeps close to the ligands but shows fluctuations at the end of the simulation. These fluctuations after adsorption indicate that the adsorption is not as stable as the one observed for HSA. The higher bHb−ligand interaction energy compared to the HSA−ligand energy (see Figure 2) underlines this. Only the interaction of O2 is strong enough for bHb to adsorb; for the other initial orientations further adjustment would be needed. From the weak interaction between bHb and the ligands and the complexity of the movements, it can be speculated that the initial orientations, which do not lead to adsorption during this simulation, may lead to adsorption after longer simulation times or bHb may finally move away from the ligands without adsorption (as for O6). The adsorption capacity of bHb is lower compared to that of HSA, and bHb rotates much stronger than HSA; for example, O3 rotates more than 250°. This is because bHb’s isoelectric point equals 7.0,40 its overall charge is low at pH 7.0, and hence there is no strong electrostatic interaction between bHb and the ion-exchange media. Rotation is needed to get a favorable orientation for the adsorption of bHb. The trend of a decreasing center of mass distance indicates that the rotation adjustment of bHb leads to a stronger interaction with the surface. Slightly charged bHb
adsorbs onto Q Sepharose FF slowly and weakly. Similar results were obtained from experiments.31 Comparing the adsorption behaviors of HSA and bHb from Figures 2−4, both HSA and bHb can adsorb onto Q Sepharose FF. Due to stronger interactions with the ligands, HSA adsorbs in preference to bHb and shows stronger adsorption with fewer fluctuations than bHb, indicated by the lower HSA−ligand interaction energy. It can also be observed from the results that the initial protein orientation influences the adsorption, which needs to be considered, especially for the weakly adsorbed proteins, and has also been discussed by previous researchers.23,24,41 Of particular interest are the residues that are close to the surface. The distance of every amino acid (AA) with the ligand surface over time is shown in Figure 5 for HSA and in Figure 6 for bHb. Which part of a protein interacts with the ligands can be extracted from these distance maps. Since the number of binding sites υ (characteristic charge) were calculated (see Table 2) as defined in section 4, the υ residues having a ligand distance of
