A molecular modeling based method to predict elution behavior and binding patches of proteins in multimodal chromatography

Accepted Manuscript Title: A Molecular Modeling Based Method to Predict Elution Behavior and Binding Patches of Proteins in Multimodal Chromatography Authors: Suvrajit Banerjee, Siddharth Parimal, Steven M. Cramer PII: DOI: Reference:

S0021-9673(17)30927-5 http://dx.doi.org/doi:10.1016/j.chroma.2017.06.059 CHROMA 358629

To appear in:

Journal of Chromatography A

Received date: Revised date: Accepted date:

26-1-2017 7-6-2017 22-6-2017

Please cite this article as: Suvrajit Banerjee, Siddharth Parimal, Steven M.Cramer, A Molecular Modeling Based Method to Predict Elution Behavior and Binding Patches of Proteins in Multimodal Chromatography, Journal of Chromatography Ahttp://dx.doi.org/10.1016/j.chroma.2017.06.059 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A Molecular Modeling Based Method to Predict Elution Behavior and Binding Patches of Proteins in Multimodal Chromatography Suvrajit Banerjee,†,§ Siddharth Parimal,†,§ Steven M. Cramer*,†,§ †

Howard P. Isermann Department of Chemical and Biological Engineering, Rensselaer Polytechnic

Institute, Troy, NY 12180, USA §

Center for Biotechnology and Interdisciplinary Studies, Rensselaer Polytechnic Institute, Troy,

NY 12180, USA *To whom correspondence should be addressed: Email: [email protected] A Molecular Modeling Based Method to Predict Elution Behavior and Binding Patches of Proteins in Multimodal Chromatography Suvrajit Banerjee, Siddharth Parimal, Steven M. Cramer Highlights 









A computationally efficient combination of short molecular dynamics simulations and continuum solvent based coarse-grained free energy calculations was employed to calculate binding free energy of proteins to surfaces presenting multimodal (MM) chromatographic ligands as a function of polar and azimuthal angles of rotation of the protein. Subsequently an iterative approach was developed to establish correlations between overall free energies of binding (average binding free energy of the most strongly bound structures) and retention factors from isocratic chromatography experiments using two model proteins, Chymotrypsinogen-A and Horse Cytochrome C on Capto MMC and Nuvia cPrime MM surfaces. These correlations along with analytical conversion relationships from the literature were used to predict protein gradient elution salt concentrations as well as selectivity reversals for two other proteins, Ubiquitin and Ribonuclease-A on the aforementioned MM surfaces. Finally, patches on protein surfaces that interacted strongly with MM surfaces were identified by determining the frequency of heavy atom contacts with the atoms of the MM SAMs. It was observed that protein faces possessing charge complementary to that of the MM surface was the main driver for binding, however, the binding was accentuated significantly by the presence of hydrophobic patches. The ability to predict elution trends and to identify key binding patches on proteins may have significant impact on process development for the separation of bioproduct related impurities.

1

Abstract Multimodal (MM) chromatography provides a powerful means to enhance the selectivity of protein separations by taking advantage of multiple weak interactions that include electrostatic, hydrophobic and van der Waals interactions.

In order to increase our

understanding of such phenomena, a computationally efficient approach was developed that combines short molecular dynamics simulations and continuum solvent based coarsegrained free energy calculations in order to study the binding of proteins to Self Assembled Monolayers (SAM) presenting MM ligands. Using this method, the free energies of protein-MM SAM binding over a range of incident orientations of the protein can be determined. The resulting free energies were then examined to identify the more “strongly bound” orientations of different proteins with two multimodal surfaces. The overall free energy of protein-MM surface binding was then determined and correlated to retention factors from isocratic chromatography. This correlation, combined with analytical expressions from the literature, was then employed to predict protein gradient elution salt concentrations as well as selectivity reversals with different MM resin systems. Patches on protein surfaces that interacted strongly with MM surfaces were also identified by determining the frequency of heavy atom contacts with the atoms of the MM SAMs. A comparison of these patches to Electrostatic Potential and hydrophobicity maps indicated that while all of these patches contained significant positive charge, only the highest frequency sites also possessed hydrophobicity. The ability to identify key binding patches on proteins may have significant impact on process development for the separation of bioproduct related impurities. Keywords: Multimodal chromatography; Molecular dynamics simulations; Coarse-grained free energy calculation; Elution prediction; Protein-surface interactions

1. INTRODUCTION Recent work has demonstrated that multimodal (MM) chromatographic resins can offer improved selectivities as compared to traditional single interaction modes such as ion exchange [1-14]. There are many forms of MM materials that are used in bioprocessing which includes mixed-beds [15], hydrophobic charge induction (HCIC) [16-18] and multimodal ligands, that is, multiple interaction modes on a single ligand [13, 14, 19-21]. 2

Various modeling efforts have been made over the past few decades to explain protein retention behavior in ion exchange systems. Initial attempts utilized simplified representations of the protein, such as a charged surface [22] or a sphere having a charge same as that of the net charge of the protein at the center of the sphere [23]. Salt dependence of protein retention in cation-exchange chromatography has been examined using the charge regulated slab model [24]. The Poisson-Boltzmann equation has been employed to efficiently predict the ionic strength dependence of retention [25]. The nonlinear PoissonBoltzmann equation, was used to predict binding equilibrium constants for protein ionexchange in terms of protein size, ζ potentials and electrolyte concentration [26]. Roth et al. [27] successfully delineated the electrostatic contribution to the enthalpy and entropy of protein adsorption in ion exchange. Dismer and Hubbuch have used molecular dynamics (MD) simulations to develop a protocol to find the preferred binding orientation for proteins on ion-exchange surfaces [28]. Previous efforts in our group employed coarse grained force fields to calculate protein-surface binding free energies at various orientations using static structures of proteins [29]. Using this approach, qualitative similarities were observed between these free energies and gradient elution salt concentrations in ion-exchange chromatography. Rigorous all-atoms MD simulations can provide the most realistic estimate of binding free energies. There are various methods for estimating free energies from these all-atoms molecular dynamics simulations [30], such as Thermodynamic Integration, Free Energy Perturbation, Umbrella Sampling or enhanced sampling methods like Replica Exchange, Metadynamics etc. [31-37]. Although accurate, the computational expense of these MD simulations is very high. Coarse-grained (continuum solvent-based) free energy calculations can circumvent this issue and provide an estimate of free energies with high computational efficiency. Solvation effects in these calculations are achieved by employing dielectric models such as Poisson-Boltzmann (PB) [38-40] and Generalized Born

(GB)

models

[41-43].

In

addition,

the

molecular

mechanics/Poisson-

Boltzmann/surface area (MM-PBSA) force field [44] which combines gas phase interaction energies with PB calculations and nonpolar force fields can be employed to yield binding free energies.

3

Electrostatic free energies are obtained by integrating the electrostatic potential field over the space around the macromolecule of interest. Various software packages are available for solving the PB equation in order to yield this electrostatic potential field. These include DelPhi [45], GRASP [46], MEAD [47], UHBD [48] and the PBEQ [49] module in CHARMM [50] which are based on the finite difference method as well as the APBS (Adaptive Poisson Boltzmann Solver) [38] which is based on the finite element method. Polar binding free energy calculation using continuum solvation methods requires accurate parameterization of atomic radii and partial charges (e.g. the PARSE force field [51]). Nonpolar solvation free energies are calculated using surface area and excluded volume models and models which include van der Waals interactions of the macromolecules with the surrounding solvent [52-54]. Hummer’s hydrophobic force field is another approach for estimating the nonpolar free energy of binding [55]. Theoretical efforts have been made to describe adsorption behavior in MM systems through isotherms by incorporating only electrostatics or both electrostatic and hydrophobic components into these isotherms [11, 56]. However, prediction of protein adsorption behavior in MM systems using molecular modeling has been limited to date. Our group has studied interactions of proteins with ligands freely dispersed in solution by employing MD simulations in order to predict preferred binding patches on protein surfaces and the chemical nature of such patches [57-59]. Quantitative Structure Activity Response (QSAR) models have also been employed to predict retention behavior in MM chromatography [60, 61]. The importance of in-silico tools for protein surface characterization of electrostatic potential [38] and surface hydrophobicity [62] has also been demonstrated in order to explain binding and elution trends in MM systems by our group [63, 64]. In addition to these studies from our group, Hirano et al. studied the interactions of the amino acid arginine with Capto MMC multimodal ligand [65], Dismer and Hubbuch predicted protein retention for ion-exchange chromatography from protein three dimensional structure [66], and Insaidoo et al. utilized molecular modeling techniques to predict retention characteristics of insulin and its variants over a range of purification challenges [67]. In this work, the main goal is to establish a predictive tool to determine gradient elution trends and preferred binding faces of various proteins in multimodal

4

chromatography systems. A framework is developed using short MD simulations and coarse-grained free energy calculations to estimate protein-surface binding free energies. This method is applied to protein-MM surface systems to yield binding free energies for a range of protein-surface binding orientations. Subsequently, an iterative method is employed to obtain overall free energy of binding by establishing a correlation between the average binding free energy of all “strongly bound” orientations with experimental isocratic retention factors.. These overall free energies show linear correlation with isocratic chromatographic retention factors. Analytical conversion relationships are then used to predict gradient elution trends for other proteins in these MM systems where predicted selectivity reversals are validated using gradient elution experiments. Finally, the strongly bound orientations are analyzed further to yield fundamental insights into protein-MM surface binding whereby we demonstrate that one particular face of each protein is mainly involved in binding and that these binding faces possess strong complementary charge to that of the MM ligand surface and the binding is further accentuated by the presence of hydrophobic patches on the preferred face.

2. THEORETICAL SECTION 2.1. Electrostatic solvation free energy ( Gelectrostatic ) The Poisson-Boltzmann (PB) equation [38, 68, 69] is solved to generate the electrostatic potential field around a macromolecule and, subsequently, this field is integrated over the space around the macromolecule to yield the electrostatic solvation free energy. The PB equation is given as follows: 4 f  r   r  r    r  r   r  sinh[ r ]  0 kT 2

(1)

Here, r denotes spatial coordinates in the volume where calculations are being done,  r  is the solvent dielectric constant,  r  is the electrostatic potential,  f  r  are the molecular charges of the protein,  r  is the ion accessibility function at r and is equal to 0 inside molecule of interest and 1 otherwise, k is the Boltzmann constant, T is absolute

5

temperature.  r 2 is the reciprocal of square of Debye length lD which, in turn, is related to the salt concentration of the solution:  r 2  8

c

 2 k k

z q2

k

2 r kT



1 . ck is the bulk 2 lD

concentration of ion k , zk is its valency and q is the electronic charge. The salt concentrations for both ions enter Equation 1 through the above expression for  r 2 , thus, enabling us to calculate  r  and, subsequently, the Gelectrostatic at different salt concentrations. In this work, the following approximation is used to convert the above nonlinear partial differential equation into a linear one (see Equation 2) under the assumption that the potential is small: sinh[ r ]   r 

(2)

Electrostatic free energy of solvation or the free energy of charging a macromolecule in solution is given by Equation 3 below. Gelectrostatic 

 

1 f r  r dr   2

(3)

2.2. Nonpolar solvation free energy ( Gnonpolar ) A model for determining the free energy change upon binding due to nonpolar interactions Gnonpolar is given by Equation 4 below [54]:

6

Gnonpolar  A  pV  Gsolvation,vdW

(4)

where, N

Gsolvation,vdW  bulk   ui( att ) ( xi , r ) x , r )dr i 1 

(5)

Here x denotes spatial coordinates of atoms of the entity under question (protein, surface or complex), r denotes spatial coordinates in the volume where calculations are being done, A is the solvent accessible surface area, p is a scaling constant accounting for mechanical work done to desolvate the cavity formed by the entity under question (protein, SAM or complex), V is the cavity volume occupied by this entity, bulk is the bulk solvent density, N is the total number of atoms of the entity under question, ui( att ) ( xi , r ) is the attractive portion of the Lennard-Jones potential of ith atom of the entity under question,  x, r ) is the Heaviside step function which is 0 for any point in the system where solvent

is excluded and 1 otherwise. The first two terms in Equation 4 are free energies of solvation due to formation of solvent excluded surface area and volume by the entity under question and the third term accounts for the attractive van der Waals interactions between the entity and surrounding solvent. 2.3. Orientation Dependent Protein-Surface Binding Free Energy [ G    ]

7

The protein-surface binding free energy Gbinding in solution (that is, the change in excess chemical potential of the protein when it is near the surface and separated from it by an infinitely large distance), is given by Equation 6 [70]. Gbinding  Gcomplex  (Gprotein  Gsurface )

(6)

Gcomplex , G protein and Gsurface are the solvation free energies of the protein-surface complex,

the protein and the surface, respectively. Each individual component (G) can be further broken down into three different components as given by the three terms in Equation 7. G  Gelectrostatic  Gnonpolar  TSconfig

(7)

Gelectrostatic and Gnonpolar refer to the electrostatic and nonpolar components of solvation free

energy, respectively, and methods to estimate these terms are delineated in Sections 2.1 and 2.2. The TSconfig term refers to the configurational entropy of the entity under question. Using Equations 3, 4, 6 and 7, the following equation for binding free energy is obtained: Gbinding  (Gelectrostatic  A  pV  Gsolvation,vdW )  T Sconfig

(8)

It has been shown that the entropy change upon binding in MM systems is primarily due to favorable gain in desolvation entropy [63] which is to be captured through the terms in parenthesis in Equation 8 as long as a sufficiently accurate bound structure can be generated using short MD simulations. Hence, here it is assumed that configurational entropy loss will have a very negligible role to play in overall binding energetics. In other words, from Equation 8 and above assumption, we conceptualize that the sum of changes in

the

different

components

of

8

solvation

free

energies

only

(Gelectrostatic  A  pV  Gsolvation,vdW ) , denoted here as simply G (see Equation 9),

should reflect the differences in binding strengths of various proteins on a given surface. G  Gelectrostatic  A  pV  Gsolvation,vdW

(9)

It is important to note that the MM-PBSA force field [44] instructs that one also include gas phase Coulombic and van der Waals interactions in Equation 9. However, we assume here that in a gas medium, Columbic interactions will be negligible as the protein and surface would exist in their neutral states as salts. The gas phase van der Waals interaction energy for spherical particles with surfaces is directly related to the size of the particles [71]. The gas phase van der Waals free energies are assumed here to be of similar values for all proteins on a given surface due to all the studied proteins being roughly globular and of similar size. Hence, such gas phase terms are not included in binding free energy calculations in the present work and the changes in solvation free energy (Equation 9) are assumed to mainly contribute to the differences in free energy of binding of various proteins on a given MM surface. G will depend on the orientation of the protein with respect to the surface. Hence, here

we consider the dependence of G on the polar and azimuthal angles of rotation of the protein     (which are described in the Methods section) as follows: G    Gelectrostatic    A     pV     Gsolvation,vdW   

(10)

2.4. The overall binding free energy GTotal and its relationship with experimental retention factors from isocratic chromatography When operating in the linear part of the protein-surface adsorption isotherm, the isocratic capacity factor k’ is defined as: k'

t R  t0  K eq t0

(11)

9

Here t R is the isocratic retention time of the protein of interest, t0 is the time taken by an unretained solute to pass through the column, K eq is the equilibrium constant of adsorption,  is the column phase ratio (i.e., surface area available for adsorption per unit mobile phase volume). Taking natural logarithm on both sides of Equation 11, we obtain, ln k '  ln Keq  ln 

(12)

The overall free energy of binding GTotal at a given salt concentration is:

GTotal   RT ln Keq

(13)

Here R is the gas constant and T is the absolute temperature. Combining Equation 12 and Equation 13, the following relationship between calculated GTotal and experimental isocratic retention factors k ' is obtained, GTotal  RT ln k ' RT ln 

(14)

10

At a particular salt concentration, for a given protein-surface pair, we assume that GTotal is the average of G    of all “strongly bound” orientations (defined further below in this section). That is,

GTotal  G   

(15)

strongly bound

Equation 14 can be rewritten as follows. GTotal  RT ln k ' intercept

(16)

Equation 16 shows that a fit of GTotal versus ln k ' should yield a straight line with a slope equal to RT and an intercept of intercept (equal to  RT ln  ). Phase ratio  is greater than unity. Hence, intercept is expected to be negative (see Equation 16). However, the exact values of the slope and the intercept are not important. This is because

GTotal  G   

strongly bound

calculated using coarse-grained approaches will not be

equal to but directly correlated with true binding free energies as has been shown by multiple authors (see reference [72-74]). Thus, according to Equation 16, on a given surface, a fit of GTotal versus ln k ' should still be linear with a positive slope. It is important to note here that the size of proteins will affect the phase ratio because larger proteins may not be able to access all of the pores in the chromatographic beads due to size exclusion effects. However, for a class of proteins of roughly the same size and a given MM surface, intercept will be constant. Since pore radius for Capto MMC is almost 10 times that of the radius of the largest protein studied here (Chymotrypsinogen A) and that of Nuvia cPrime even larger [60, 75], all pores are assumed to be equally accessible by the proteins studied resulting in a constant phase ratio. We assume that for a given MM surface, the protein orientations which exhibit G    < 0 kJ/mol at a particular low salt concentration Clow are the ones primarily responsible for strong protein-surface binding and, hence, are the “strongly bound” orientations mentioned

11

above. Further, we assume that these orientations will play the most important role in dictating protein elution behavior at higher salt concentrations. Thus, only these structures are included in free energy averaging for all higher salt concentrations using Equation 15. The method for determining Clow for a given MM surface is described in the Materials and Methods section. 2.5. Prediction of gradient elution salt concentration The stoichiometric displacement model (SDM) is commonly used for correlating isocratic retention factors with the ionic strength of solution: ln k '  Z ln I 1  ln K

(17)

Z is the number of binding sites, K is a constant and I is the ionic strength [76, 77].

Since, in this case, a 1:1 electrolyte (NaCl) is used, the ionic strength is equal to salt concentration. According to Equation 17, a linear fit of ln k ' versus ln I 1 will yield the quantities K and Z which can then be employed to predict the gradient elution salt concentration by using the relationship developed by Parente and Wetlaufer [78]. For a gradient with initial starting salt concentration of 0M, this relationship is as follows:  1   

I Gradient

V K ( Z  1) I Final   Z 1   0  VG  

(18)

I Gradient is the gradient elution salt concentration for the protein under question, V0 is the

dead volume for an unretained solute to elute, I Final is the final gradient salt concentration and VG is the total volume of the gradient. 3. MATERIALS AND METHODS 3.1. In-silico representation of proteins and surfaces Model Self Assembled Monolayer (SAM) surfaces presenting MM ligands were constructed. The SAMs had one C10 strand attached to a base sulfur atom. The sulfur atoms were position-restrained to locations corresponding to those on gold 111 lattice [79]. The alkane tails were represented using the OPLS united atom representation [80]. MM ligand 12

head groups used for this study are Capto MMC and Nuvia cPrime which are depicted in Figure 1. Parameters for these ligands were assigned using parameters of constituent groups from the all-atom Assisted Model Building with Energy Refinement (AMBER) Parm-94 force field [57, 81]. This force field has been employed previously in our research group to parameterize ubiquitin and Capto MMC ligands for free ligand simulations and proved to be effective in predicting the binding hotspots of these ligands on ubiquitin that were observed in analogous NMR experiments [58]. Using these head groups, a SAM surface comprising Capto MMC and hydroxyl head groups on each SAM slab (surface density of Capto head groups was 1.14 ligands per nm2). This surface coverage density closely resembled the chromatographic surface ligand density [82]. Water molecules near –OH surfaces were previously shown to exhibit behavior akin to that of bulk water [83] and proteins also exhibit weak binding to such surfaces [84]. Hence, these hydrophilic hydroxyl head groups offered a weakly interacting background to enable interactions of the protein with the MM ligands only. The proteins ubiquitin (pdb ID: 1D3Z), chymotrypsinogen A (2CGA), ribonuclease A (1RBX) and horse cytochrome C (1HRC) were modeled using Amber94 force-field [81]. Charges were assigned using PROPKA [85] at pH 5.0. Suitable number of Na+ counterions was added to neutralize the system. Water molecules were represented using the extended TIP3P model. 3.2. Generation of bound structures of proteins to MM SAM surfaces at various sets of incident angles using MD simulations To obtain overall free energy of binding it is first necessary to generate bound structures of proteins to MM surfaces. For this purpose, 2 nanosecond long MD simulations were performed over a set of incident angles of the protein,  and  (See Figure 2). These sets of angles     represented points equidistant across the surface of a unit sphere and 101 such points were chosen. The starting structure for each orientation was chosen such that the separation distance between the lowest protein atom and the highest SAM atom was 0.45 nm. To maintain a specific orientation of the protein to the SAM, the C- atoms of three core residues of each protein were position restrained in the plane parallel to plane of the surface. The main idea behind performing these short MD simulations was to obtain a reasonable estimate of the bound structure of the protein to the SAM at each orientation. 13

All-atoms molecular dynamics simulations were performed using GROMACS [86-88]. All simulations were conducted in the isothermal-isobaric ensemble at a temperature of 298K and a pressure of 1 atm using the Nose-Hoover thermostat [89, 90] (time constant for coupling was 0.5 ps). An anisotropic Parrinello-Rahman barostat [91] (time constant was 1.0 ps) was employed to allow independent relaxation only in the direction perpendicular to the SAM surface. The Particle Mesh Ewald method [92] was used to calculate electrostatic interactions. Real space cutoff and cutoff for Lennard-Jones interactions were both taken to be 1 nm. A time step of 2 fs was used. Configurations were saved every 1 ps for analysis. By performing these simulations, proteins bound to the surface if a favorably interacting face was exposed to the MM SAM and were repelled from the surface if a face was exposed to the SAM that had unfavorable interactions (exemplified in Figure S1A and Figure S1B). After these bound structures were obtained, the water and ions were removed and continuum based coarse-grained binding free energy calculations were performed as described below. The free energies obtained for one orientation of ubiquitin and horse cytochrome C on Capto MMC for two different salt concentrations as a function of simulation time is shown in Figure S1C. As can be seen in this Figure, the free energies remained largely invariant after 2 ns of simulation, showing that this simulation time was sufficiently long for obtaining a stable bound structure. 3.3. Continuum solvent based coarse-grained free energy calculations Prior to performing coarse grained free energy calculations, water molecules and ions were removed from around the protein-MM SAM bound structures generated as mentioned in the previous section. The protein and/or SAM were then re-parameterized using partial charges and atomic radii prescribed in the PARSE force-field [51] that are continuum solvent-based. The linearized Poisson-Boltzmann equation (that is, Equation 1 with the modification depicted in Equation 2) was solved using the Adaptive Poisson-Boltzmann Solver software (APBS) package [38] on a coarse grid having a grid thickness of ~0.08 nm in each direction centered on the molecule of interest (that is, protein, MM SAM or proteinSAM complex). This calculation was further refined through a secondary calculation on a smaller grid focusing on the protein-SAM binding region with a grid thickness of ~0.03 nm. The solvent and protein/SAM dielectric constants were taken to be 78.54 and 2.0,

14

respectively. The mobile ion charges were taken to be +1 and -1 and atomic radii 0.095 nm and 0.181 nm [69] to emulate sodium and chloride ions, respectively. Electrostatic solvation free energy was calculated using Equation 3. The nonpolar solvation free energy was calculated using Equation 4. The surface areas and volumes of protein, surface and complex for nonpolar free energy calculations were calculated using the g_sas module of GROMACS [88]. The van der Waals solvation free energy was calculated using an inhouse code. The values of the various parameters employed in calculating the nonpolar free energy of binding can be found in reference [54]. G    ’s were then determined for various sets of     using Equation 10. Here, it was assumed that the bound structures will not alter significantly such that the structures generated at low salt using MD simulations can be used for evaluating G    at higher salt concentrations. 3.4. Overall free energies of binding and its correlation with chromatography experiments In this work, an iterative scheme was developed for calculating the overall free energy of binding ( GTotal ) by establishing a linear correlation between the calculated GTotal ’s and isocratic retention factors ( k ' ) for a given MM surface. This iterative process was developed to determine the Clow concentration which in turn determined the “strongly bound” orientations (see Theory section). Starting from an initial value of Clow =0.1M, the “strongly bound” structures (that is, bound structures which exhibit G    < 0 kJ/mol) were determined for two proteins, chymotrypsinogen-A and horse cytochrome C, on a given MM surface. The GTotal ’s were then calculated using Equation 15 for both proteins at four higher salt concentrations. A plot of the calculated GTotal versus ln k ' for these two proteins at the four higher salt concentrations was subsequently generated and the R 2 of the corresponding linear fit was determined (As discussed in the Theory section, Equation 16 depicts that a plot of GTotal versus ln k ' should be linear with a positive slope). This process was repeated at successively higher Clow ’s until a final Clow value was obtained which yielded the minimum error (i.e., maximum R2) in the linear correlation.

15

This final value was regarded as the Clow value for the given MM surface, and the linear fit of GTotal versus ln k ' corresponding to this Clow was established as the correlation to be used for converting calculated GTotal ’s to theoretical isocratic k’s for other proteins on that MM surface. 3.5. Prediction of gradient elution salt concentration Once the correlation between GTotal and ln k ' was established for a given MM surface, GTotal ’s were determined for two other proteins, ubiquitin and ribonuclease A, at different

salt concentrations by averaging G    ’s for all “strongly bound” orientations (Equation 15) determined at the Clow for that MM surface, and the correlation was employed to predict theoretical isocratic k ' s. Next, a fit of these theoretical k ' s versus corresponding salt concentrations on a logarithmic plot yielded constants K and Z for a protein binding to a given surface (see Equation 17). Finally, constants K and Z coupled with chromatographic column parameters were substituted into Equation 18 to predict gradient elution salt concentrations for a protein on a given MM system. 3.6. Identifying preferred binding regions on the protein surface The “strongly bound” orientations obtained from the Clow calculations (described in the previous section) were analyzed further to determine patches on the protein surface that were involved in binding. For this purpose, the frequency of protein heavy atoms in contact with MM SAM atoms was determined from all the “strongly bound” structures. A protein heavy atom located within 0.6 nm of any atom of the MM SAM was considered to be in contact with the SAM (~0.7 nm is the solvent separated methane-methane distance in water [93]), and was assigned a count of 1 for every instant it fulfilled this criterion. Finally, the counts were normalized and mapped onto the protein surface (color scheme discussed in Results and Discussions). Electrostatic potential maps [38] and Spatial Aggregation Propensity maps [62] (hydrophobicity maps) were also generated to examine qualitatively the chemical nature of the strongly binding patches. 3.7. Isocratic and gradient chromatography experiments

16

Capto MMC was purchased from GE Healthcare (Uppsala, Sweden) and Nuvia cPrime was obtained from Bio-Rad Laboratories (Hercules CA, USA). These chromatographic media were packed into Pharmacia Biotech 1ml glass columns. Hydrochloric acid, sodium chloride, acetic acid and sodium acetate, horse cytohrome C, chymotrypsinogen A and ribonuclease A were purchased from Sigma-Aldrich (St. Louis, MO). Sodium hydroxide was purchased from Thermo Fisher Scientific (Pittsburgh PA, USA). Ubiquitin was purchased from Boston Biochem (Cambridge MA, USA). Chromatography experiments were performed on an Äkta Explorer 100 (Amersham Biosciences, Uppsala, Sweden). Multimodal columns were equilibrated with 30 column volumes of 20 mM sodium acetate (pH 5.0) prior to the start of every experiment and regenerated with 5 column volumes of 1M sodium hydroxide. In the isocratic chromatography experiments, salt of desired concentration was also present in the equilibration buffer. Isocratic

chromatography

experiments

were

performed

using

proteins

chymotrypsinogen A and horse cytochrome C. These experiments were carried out with constant mobile phase NaCl compositions (0.6M, 0.65M, 0.7M and 0.75M for Capto MMC, and at 0.3M, 0.35M, 0.4M and 0.45M for Nuvia cPrime). All experiments were performed at a flowrate of 1 mL/min. The buffer consisted of sodium acetate (20 mM sodium acetate) adjusted to pH 5.0. The column effluent was monitored at a UV wavelength of 280 nm. Retention times t R were determined from the peak maximum of each protein sample. Linear gradient chromatography experiments were carried out with the proteins ubiquitin and ribonuclease A using a linear salt gradient from buffer A (20 mM sodium acetate, pH 5.0) to buffer B (20 mM sodium acetate, 1.5 M sodium chloride, pH 5.0) over 45 column volumes at a flowrate of 1 mL/min. Injection volume was 100 μL. The column effluent was monitored at a UV wavelength of 280 nm. Retention times were determined from the peak maximum of each protein sample. Elution salt concentrations were calculated from the retention time data using the specified linear salt gradient conditions. Chromatographic parameters corresponding to Equation 18 (above) for these columns were: IFinal=1.5M for both columns and

V0 = 0.029 and 0.023 for Capto MMC and Nuvia VG

cPrime, respectively.

17

4. RESULTS AND DISCUSSIONS 4.1. Protein-surface interaction maps Binding free energy maps were generated using the approach described above for a range of incident angles [various sets of  ,   ] for different proteins on a given multimodal surface. G  ,   at each point on these contour plots were calculated from protein-MMC SAM bound structures using Equation 10 after performing 2 ns of MD simulation. For all of the MD simulations, the protein faces that bound to the MM surface were seen to only interact with MM head groups and maintained an average distance between the hydrophilic –OH background of the SAM and the closest protein atom of ~1 nm. Representative binding maps are shown in Figure 3 for the binding of the 4 proteins evaluated in this study (-chymotrypsinogen A, horse cytochrome C, ubiquitin and ribonuclease A) to the Capto MMC SAM surface in 0.25M NaCl at pH 5.0. The binding free energy scale is depicted up to -20 kJ/mol which corresponds to quite strong interactions (8kT). The binding maps in Figure 3 indicate that there were a relatively few number of highly preferred orientations as compared to unfavorable ones. As will be discussed later, these favorable orientations generally tend to be associated with a single face of the protein surface. Figure 3 also depicts the protein-SAM bound structure for the strongest protein binding orientation corresponding to each map. 4.2. Correlation between calculated overall binding free energies GTotal and experimental retention factors from isocratic chromatography As described in the Methods section, an iterative process was developed to establish a linear correlation between GTotal values and the experimental isocratic retention factors (k’) on a given multimodal resin system (see Section 3.4). This approach determined an appropriate low salt concentration ( Clow ) at which the protein orientations that contributed to strong binding and their associated bound structures on a given surface could be selected for optimal prediction of GTotal at higher salt concentrations. As mentioned in the Theory and Methods sections, the assumption here is that for a given protein and resin combination, only these “strongly bound” protein orientations will affect retention at

18

higher salt concentrations (see Theory Section 2.4 and Methods Section 3.4 for more details). The proteins employed in this study for developing the correlations on the two multimodal resin systems were chymotrypsinogen-A and horse cytochrome C which possess different levels of surface hydrophobicities and charge characteristics. Using this iterative method, it was determined that for Capto MMC, a Clow of 0.25 M NaCl resulted in the highest R2 (R2=0.91) value for a linear correlation between GTotal and ln k ' over a range of higher salt concentrations (see Figure 4A). The slope of the linear

correlation plot was 1.55 kJ/mol and the intercept was -1.17 kJ/mol. A similar protocol was followed to find the linear relationship between calculated GTotal with experimental isocratic k’ for these two proteins on the Nuvia cPrime (R2=0.93) surface and Figure 4B shows the resulting linear fit. The Clow value obtained here was found to be 0.14 M. This is not surprising considering that proteins are generally less retained on Nuvia cPrime as compared to Capto MMC [60]. Interestingly, in previous work from our group with these multimodal resins, it was found that all of the basic proteins studied eluted above these particular salt values for the two resins. The slope of the linear correlation plot for Nuvia cPrime was 1.45 kJ/mol and the intercept was 0.37 kJ/mol. To summarize, the linear correlations between GTotal and ln k ' were as follows: For Capto MMC, GTotal  1.55ln k '1.17

(19a)

For Nuvia cPrime, GTotal  1.45ln k ' 0.37

(19b)

According to Equation 16, the slope of these linear fits was expected to be equal to RT (that is, equal to 2.48 kJ/mol at 298K). Although the slopes were less, it is encouraging that, not only do they lie in the same order of magnitude as that of the theoretical slope, but also differed from the theoretical value by < 1 kJ/mol. Figure 5 illustrates how the

G  ,   of the “strongly bound” orientations for these two proteins on Capto MMC change with increasing salt concentration. (Note: Each map in Figure 5 corresponds to a single GTotal value on the linear correlation plot in Figure 4A). 4.3. Prediction of gradient elution salt concentration

19

Once the Clow value and the relationship between calculated binding free energies GTotal and experimental isocratic capacity factors k’ was established for the two resins, a theoretical framework was employed to directly predict gradient elution trends for two additional proteins (ubiquitin and ribonuclease A) without any further experimental data (see Theory Section 2.5 and Methods Section 3.5). The determination of Clow enabled the identification of the “strongly bound” orientations, which in turn was used to determine the binding maps at higher salt concentrations as shown in Figure

S2. At each salt

concentration, the GTotal was then determined using Equation 15 and the previously developed linear relationship for each resin (Equations 19a and 19b) was employed to predict the isocratic k’ values at each salt concentration. Equation 17 was then used to estimate parameters K and Z for these proteins binding to the Capto and Nuvia surfaces (linear fits are shown in Figure S3 and the fit parameters are listed in Table S1). It is encouraging that the data in Figure S3 which is generated using the present approach was linearly very well correlated (R2 is ~0.99 for all fits in this Figure). Finally, these fit parameters K and Z were used in Equation 18 to predict gradient elution salt concentrations for the two proteins on the two MM ligand surfaces. The predicted and experimental elution salt concentrations are depicted in Figure 6. As can be seen in Figure 6, the gradient elution salt concentrations were predicted within 30mM of experimental values for the Capto MMC and within 100 mM for the Nuvia cPrime resin systems. Importantly, this predictive approach was able to capture the selectivity reversal of these two proteins on the Capto MMC and Nuvia cPrime resin surfaces. This is an important achievement of the developed method. The pI of ubiquitin is 6.8 and that of ribonuclease A is 9.6. At the pH employed here (pH 5.0), one might expect earlier elution of ubiquitin than ribonuclease A on these cation-exchange systems. However, Capto MMC is structurally very different from Nuvia cPrime (see Figure 1) and hence can interact with proteins via mechanisms different from that of Nuvia cPrime. The fact that our method captures this qualitatively different elution behavior is highly encouraging and indicates that this combination of short MD simulations and coarse grained free energy calculations has the potential to delineate unique selectivity trends offered by different multimodal chromatographic systems.

20

4.4. Preferred binding regions on protein surface Not only does the approach presented in this paper help to determine the binding energetics of various protein orientations to the multimodal surface, but it also enables one to determine the relative importance of various protein surface patches involved in the binding. Using the method described in Section 3.6, the frequency of various protein heavy atoms binding to the multimodal SAM for the “strongly bound” structures was calculated and the proteins were colored appropriately. For each protein, the colored structures are presented for two opposite faces of the protein associated with binding to the Capto MMC and the Nuvia cPrime surfaces. In addition, the electrostatic potential (EP) maps and the spatial aggregation propensity (SAP) maps are presented for the two protein faces. Finally, predominant patches on the protein surfaces resulting from binding to Capto MMC only are also circled on the other structures to facilitate comparison. For chymotrypsinogen-A, a strong binding patch (A1) was observed on the front face for binding to both MM surfaces (Figure 7A and B). As can be seen in Figure 7C and D, this region possessed both moderate positive charge as well as hydrophobicity. The other binding patch for chymotrypsinogen-A (A2) on Capto MMC covered a wider region on the protein surface, but with lower frequencies than A1. Interestingly, this patch established more protein-MM surface heavy atom contacts with Capto MMC (see Figure 7A) as compared to Nuvia cPrime (Figure 7B). As can be seen in Figure 7C and D, the A2 region while having significant positive charge had much less hydrophobicity than A1. Even though there is a significant hydrophobic patch on the back face of this protein (see Figure 7D), this face is predominantly negatively charged (Figure 7C) and, thus the protein would not be expected to bind through this face. It is encouraging that our simulation method captures this behavior with no binding patches on the back face (Figure 7A and B). For horse cytochrome C, a strong binding patch was observed with both resins at C2 (Figures 8A and B). As can be seen in Figures 8C and D, this patch was associated with both a positively charged and significantly hydrophobic region. An additional binding patch was also observed for Capto MMC (C1) (and to a less extent with Nuvia cPrime) which was associated with a large, primarily positively charged region (Figure 8C). The back face of this protein was predominantly negatively charged (Figure 8C) and, thus the 21

protein would not be expected to bind through this face. Again, our simulation method was able to capture this behavior with minimal binding patches on the back face (Figure 9A and B). For ribonuclease-A, two binding patches were obtained for the front face of the protein (R1 and R2) on both MM surfaces. These patches were associated with positively charged regions (Figure 9C) and minimal hydrophobicity (Figure 9D). As was observed for the previously discussed proteins, the Capto MMC surface showed higher frequencies of heavy atom contacts with the protein than was observed with Nuvia cPrime. While the back face of the protein had some positively charged regions, it possessed predominantly negative electrostatic potential. Further, even though this back face of the protein was slightly more hydrophobic than the front face (Figure 9D), these hydrophobic patches tended to overlap primarily with negatively charged patches (Figure 9C). For these two reasons this protein would not be expected to bind to the MM surfaces through this face, which was corroborated using the simulation method (Figure 9A and B). For ubiquitin, two binding patches (U1 and U2) were observed for the Capto MMC surface (see Figure 10A). These patches were associated with positively charged regions and varying levels of hydrophobicities (Figure 10C and D). Binding patch U1, which encompassed both the front and back faces of the protein, was predominantly positively charged (see Figure 10C) and moderately hydrophobic in nature (Figure 10D). Patch U2 was also predominantly positively charged but possessed significant hydrophobicity (Figures 10C and D). It is interesting that the multimodal ligand with more solvent exposed hydrophobic groups (Capto MMC) (see Figure 10A) bound to this patch while Nuvia cPrime did not (Figure 10B). The U2 binding region has been previously identified in our group as an important binding face of ubiquitin to nanoparticles coated with Capto MMC ligands using Nuclear Magnetic Resonance Spectroscopy [63]. However, in that study the ligand density employed on the nanoparticle surface was three times higher than that used on the chromatographic resin system and also used in the simulations in this work, which may be the reason why U2 was not identified as being important for the Nuvia cPrime system. The importance of U1 in the current study may also be, to some extent, due to the higher flexibility of this region of the protein which could have resulted in elevated heavy atom frequency data using the current approach.

22

4.5. Computational Efficiency For a protein-surface pair, and for any one of the test proteins (that is, ubiquitin or ribonuclease-A), the current approach required 101 short MD simulations followed by Poisson-Boltzmann calculations at 5 different salt concentrations for each orientation. MD simulation of each orientation required about 20 hours to run on 64 eight-core 2.6 GHz Intel Xeon E5-2650 processors. Poisson Boltzmann calculations were performed in about ~9 hours. The nonpolar binding free energy calculations required much less time to run than Poisson-Boltzmann calculations and were also carried out in parallel with the PoissonBoltzmann calculations. Hence nonpolar free energy calculations did not limit computational efficiency. Subsequent retention predictions required minimal time. While this entire process could be achieved in approximately 30 hours, free energy calculation methods based on all-atoms molecular dynamics simulations, such as umbrella sampling, with the same protein-surface pair required significantly more time. Using the same computer hardware indicated above, all-atoms MD methods would require computational time on the order of one to two weeks to perform these calculations for just one orientation at a given salt concentration. This indicates that the current method is significantly more computationally efficient than established all-atoms MD approaches. 4.3. CONCLUSIONS In this work, a combination of short MD simulations and coarse-grained free energy calculations was used as a computationally efficient framework to determine free energy of binding of various proteins faces on MM SAM surfaces. The resulting binding maps were then examined to identify the more strongly binding orientations of the different proteins with the two multimodal surfaces. A method was also developed for calculating the overall free energy of binding GTotal at different salt concentrations for each proteinsurface pair. This method employed an iterative scheme that established a linear correlation between the calculated GTotal ’s and retention factors ( k ' ) from isocratic chromatographic experiments over a range salt concentrations. Once this correlation was established for a given resin system, GTotal ’s were evaluated for other proteins and the established correlation was used to obtain theoretical isocratic k ' s. Finally, gradient elution salt

23

concentrations were predicted using analytical expressions from these theoretical k’s. This predictive method was not only able to predict the gradient elution salt concentrations but to also capture selectivity reversals observed with Capto MMC and Nuvia cPrime for a given protein pair. Orientational preferences of proteins binding to these MM surfaces were further studied by determining the frequency of heavy atom contacts with the atoms of the SAM and coloring the protein surface appropriately. This resulted in the identification of preferred binding patches on the proteins for interacting with the two multimodal surfaces. A comparison of these patches to EP and SAP maps indicated that while all of these patches contained significant positive charge, only the highest frequency sites possessed hydrophobicity. This approach enables us to determine the relative importance of various protein orientations towards a multimodal surface with respect to their contribution to binding. It also enables a quantitative prediction of gradient elution behavior and selectivity reversals on different resin systems for a given protein pair. Further, it enables us to identify important binding patches on the protein surfaces for a given multimodal system. The work to date has focused on small proteins, however, it is important to extend this approach to larger biological products such as Fabs and mAbs. This extension will require a deeper understanding of preferred binding regions and steric limitations in these larger solute systems. This work may have significant impact on process development for the separation of bioproduct related impurities. The ability to predict whether a given MM surface would be capable of interrogating key faces of a protein where variations may occur, could be extremely useful as an initial resin screening technique and for informing the design of bioproducts for their manufacturability. ACKNOWLEDGMENTS The authors acknowledge financial support by the National Science Foundation (CBET 1160039). The authors thank Prof. Shekhar Garde for useful discussions and suggestions related to molecular simulations, statistical mechanics and adsorption thermodynamics, and the Center for Computational Innovations at Rensselaer Polytechnic Institute for computational infrastructure support. The authors also thank Julie Robinson, Swarnim Ranjan and Ronak Gudhka for experimental support.

24

REFERENCES 1.

W.R. Melander, Z. El Rassi, C. Horvath, Interplay of hydrophobic and electrostatic interactions in biopolymer chromatography : Effect of salts on the retention of proteins. J. Chromatogr. A 469 (1989) 3-27.

2.

W.S. Hancock, J.T. Sparrow, Use of Mixed-Mode, High-Performance LiquidChromatography for the Separation of Peptide and Protein Mixtures. J. Chromatogr. 206 (1981) 71-82.

3.

L.W. Mclaughlin, Mixed-Mode Chromatography of Nucleic-Acids. Chem. Rev. 89 (1989) 309-319.

4.

B.L. Johansson, M. Belew, S. Eriksson, G. Glad, O. Lind, J.L. Maloisel, N. Norrman, Preparation and characterization of prototypes for multi-modal separation aimed for capture of positively charged biomolecules at high-salt conditions. J. Chromatogr. A 1016 (2003) 35-49.

5.

B.L. Johansson, M. Belew, S. Eriksson, G. Glad, O. Lind, J.L. Maloisel, N. Norrman, Preparation and characterization of prototypes for multi-modal separation media aimed for capture of negatively charged biomolecules at high salt conditions. J. Chromatogr. A 1016 (2003) 21-33.

6.

D. Gao, D.Q. Lin, S.J. Yao, Protein adsorption kinetics of mixed-mode adsorbent with benzylamine as functional ligand. Chem. Eng. Sci. 61 (2006) 7260-7268.

7.

S.C. Burton, N.W. Haggarty, D.R.K. Harding, One step purification of chymosin by mixed mode chromatography. Biotechnol. Bioeng. 56 (1997) 45-55.

8.

S.C. Burton, D.R.K. Harding, High-density ligand attachment to brominated allyl matrices and application to mixed mode chromatography of chymosin. J. Chromatogr. A 775 (1997) 39-50.

9.

X.D. Liu, C. Pohl, A new weak anion-exchange/reversed-phase mixed-mode stationary phase for simultaneous separation of basic, acidic and neutral pharmaceuticals. LC GC Europe (2007) 33-33.

10.

K. Kallberg, K. Becker, L. Bülow, Application of a pH responsive multimodal hydrophobic interaction chromatography medium for the analysis of glycosylated proteins. J. Chromatogr. A 1218 (2011) 678-683.

25

11.

K. Zhao, L. Yang, X. Wang, Q. Bai, F. Yang, F. Wang, Preparation of a novel dualfunction strong cation exchange/hydrophobic interaction chromatography stationary phase for protein separation. Talanta 98 (2012) 86-94.

12.

K.A. Kaleas, C.H. Schmelzer, S.A. Pizarro, Industrial case study: Evaluation of a mixed-mode resin for selective capture of a human growth factor recombinantly expressed in E. coli. J. Chromatogr. A 1217 (2010) 235-242.

13.

B.L. Johansson, M. Belew, S. Eriksson, G. Glad, O. Lind, J.L. Maloisel, N. Norrman, Preparation and characterization of prototypes for multi-modal separation media aimed for capture of negatively charged biomolecules at high salt conditions. J. Chromatogr. A 1016 (2003) 21-33.

14.

B.L. Johansson, M. Belew, S. Eriksson, G. Glad, O. Lind, J.L. Maloisel, N. Norrman, Preparation and characterization of prototypes for multi-modal separation aimed for capture of positively charged biomolecules at high-salt conditions. J. Chromatogr. A 1016 (2003) 35-49.

15.

Z. el Rassi, C. Horvath, Tandem columns and mixed-bed columns in highperformance liquid chromatography of proteins. J. Chromatogr. 359 (1986) 255264.

16.

E.

Boschetti,

Antibody

separation

by

hydrophobic

charge

induction

chromatography. Trends Biotechnol. 20 (2002) 333-337. 17.

S.C. Burton, D.R. Harding, Hydrophobic charge induction chromatography: salt independent protein adsorption and facile elution with aqueous buffers. J. Chromatogr. A 814 (1998) 71-81.

18.

G. Zhao, L. Zhang, S. Bai, Y. Sun, Analysis of hydrophobic charge induction displacement chromatography by visualization with confocal laser scanning microscopy. Sep. Purif. Technol. 82 (2011) 138-147.

19.

S. Wongyai, Synthesis and characterization of phenylpropanolamine bonded silica for multimode liquid chromatography of small molecules. Chromatographia 38 (1994) 485-490.

20.

M. Lämmerhofer, R. Nogueira, W. Lindner, Multi-modal applicability of a reversed-phase/weak-anion exchange material in reversed-phase, anion-exchange,

26

ion-exclusion,

hydrophilic

interaction

and

hydrophobic

interaction

chromatography modes. Anal. Bioanal. Chem. 400 (2011) 2517-2530. 21.

O. Pitiot, L. Folley, M.A. Vijayalakshmi, Protein adsorption on histidylaminohexyl-Sepharose 4B. I. Study of the mechanistic aspects of adsorption for the separation of human serum albumin from its non-enzymatic glycated isoforms (advanced glycosylated end products). J. Chromatogr. B Biomed. Sci. Appl. 758 (2001) 163-172.

22.

J. Staahlberg, B. Joensson, C. Horvath, Theory for electrostatic interaction chromatography of proteins. Anal. Chem. 63 (1991) 1867-1874.

23.

C.M. Roth, A.M. Lenhoff, Electrostatic and van der Waals contributions to protein adsorption: computation of equilibrium constants. Langmuir 9 (1993) 962-972.

24.

E. Hallgren, F. Kálmán, D. Farnan, C. Horváth, J. Ståhlberg, Protein retention in ion-exchange chromatography: effect of net charge and charge distribution. J. Chromatogr. A 877 (2000) 13-24.

25.

B. Jönsson, J. Ståhlberg, The electrostatic interaction between a charged sphere and an oppositely charged planar surface and its application to protein adsorption. Colloids Surf., B 14 (1999) 67-75.

26.

W.R. Bowen, A.O. Sharif, Long-range electrostatic attraction between like-charge spheres in a charged pore. Nature 393 (1998) 663-665.

27.

C.M. Roth, J.E. Sader, A.M. Lenhoff, Electrostatic contribution to the energy and entropy of protein adsorption. J. Colloid Interface Sci. 203 (1998) 218-221.

28.

F. Dismer, J. Hubbuch, In silico prediction of protein binding using molecular dynamic simulations: A future tool for accelerated process development. J. Biosci. Bioeng. 108 (2009) S60.

29.

A.S. Freed, S.M. Cramer, Protein− Surface Interaction Maps for Ion-Exchange Chromatography. Langmuir 27 (2011) 3561-3568.

30.

C. Chipot, A. Pohorille (2007) Free energy calculations (Springer).

31.

S.A. Adcock, J.A. McCammon, Molecular dynamics: survey of methods for simulating the activity of proteins. Chem. Rev. 106 (2006) 1589-1615.

32.

T. Haliloglu, I. Bahar, B. Erman, Gaussiam dynamics of folded proteins. Phys. Rev. Lett. 79 (1997) 3090-3093.

27

33.

B. Isralewitz, M. Gao, K. Schulten, Steered molecular dynamics and mechanical functions of proteins. Curr. Opin. Struct. Biol. 11 (2001) 224-230.

34.

J. Ma, M. Karplus, Ligand-induced conformational changes in ras p21: a normal mode and energy minimization analysis. J. Mol. Biol. 274 (1997) 114-131.

35.

A.D. Mackerell, Jr., Empirical force fields for biological macromolecules: overview and issues. J. Comput. Chem. 25 (2004) 1584-1604.

36.

B. Raux, The Calculation of the Potential of Mean Force Using ComputerSimulations. Comput. Phys. Commun. 91 (1995) 275-282.

37.

Y. Sugita, Y. Okamoto, Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 314 (1999) 141-151.

38.

N.A. Baker, D. Sept, S. Joseph, M.J. Holst, J.A. McCammon, Electrostatics of nanosystems: application to microtubules and the ribosome. Proc. Natl. Acad. Sci. U.S.A. 98 (2001) 10037-10041.

39.

R. Luo, M.S. Head, J. Moult, M.K. Gilson, p K a shifts in small molecules and HIV protease: Electrostatics and conformation. J. Am. Chem. Soc. 120 (1998) 61386146.

40.

M. Feig, A. Onufriev, M.S. Lee, W. Im, D.A. Case, C.L. Brooks, Performance comparison of generalized born and Poisson methods in the calculation of electrostatic solvation energies for protein structures. J. Comput. Chem. 25 (2004) 265-284.

41.

A. Onufriev, D. Bashford, D.A. Case, Modification of the generalized Born model suitable for macromolecules. J. Phys. Chem. B 104 (2000) 3712-3720.

42.

W. Im, M.S. Lee, C.L. Brooks, Generalized born model with a simple smoothing function. J. Comput. Chem. 24 (2003) 1691-1702.

43.

D. Bashford, D.A. Case, Generalized born models of macromolecular solvation effects. Annu. Rev. Phys. Chem. 51 (2000) 129-152.

44.

B. Kuhn, P. Gerber, T. Schulz-Gasch, M. Stahl, Validation and use of the MMPBSA approach for drug discovery. J. Med. Chem. 48 (2005) 4040-4048.

45.

I. Klapper, R. Hagstrom, R. Fine, K. Sharp, B. Honig, Focusing of electric fields in the active site of Cu‐Zn superoxide dismutase: Effects of ionic strength and amino‐ acid modification. Proteins: Struct., Funct., Bioinf. 1 (1986) 47-59.

28

46.

A. Nicholls, Graphical representation and analysis of surface properties. Surfaces 5 (1993) 2.4.

47.

D. Bashford (1997) An object-oriented programming suite for electrostatic effects in biological molecules An experience report on the MEAD project. Scientific computing in object-oriented parallel environments, (Springer), pp 233-240.

48.

J.D. Madura, J.M. Briggs, R.C. Wade, M.E. Davis, B.A. Luty, A. Ilin, J. Antosiewicz, M.K. Gilson, B. Bagheri, L.R. Scott, Electrostatics and diffusion of molecules in solution: simulations with the University of Houston Brownian Dynamics program. Comput. Phys. Commun. 91 (1995) 57-95.

49.

W. Im, D. Beglov, B. Roux, Continuum solvation model: computation of electrostatic forces from numerical solutions to the Poisson-Boltzmann equation. Comput. Phys. Commun. 111 (1998) 59-75.

50.

B.R. Brooks, R.E. Bruccoleri, B.D. Olafson, D.J. States, S. Swaminathan, M. Karplus, CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 4 (1983) 187-217.

51.

C. Tan, L. Yang, R. Luo, How well does Poisson-Boltzmann implicit solvent agree with explicit solvent? A quantitative analysis. J. Phys. Chem. B 110 (2006) 1868018687.

52.

W. Wang, P.A. Kollman, Free energy calculations on dimer stability of the HIV protease using molecular dynamics and a continuum solvent model. J. Mol. Biol. 303 (2000) 567-582.

53.

C. Tan, Y.-H. Tan, R. Luo, Implicit nonpolar solvent models. J. Phys.Chem. B 111 (2007) 12263-12274.

54.

J.A. Wagoner, N.A. Baker, Assessing implicit models for nonpolar mean solvation forces: the importance of dispersion and volume terms. Proc. Natl. Acad. Sci. U.S.A. 103 (2006) 8331-8336.

55.

G. Hummer, Hydrophobic force field as a molecular alternative to surface-area models. J. Am. Chem. Soc. 121 (1999) 6299-6305.

56.

H.S. Karkov, L. Sejergaard, S.M. Cramer, Methods development in multimodal chromatography with mobile phase modifiers using the steric mass action model. J. Chromatogr. A 1318 (2013) 149-155.

29

57.

S. Parimal, S. Garde, S.M. Cramer, Interactions of Multimodal Ligands with Proteins: Insights into Selectivity Using Molecular Dynamics Simulations. Langmuir 31 (2015) 7512-7523.

58.

A.S. Freed, S. Garde, S.M. Cramer, Molecular simulations of multimodal ligand– protein binding: Elucidation of binding sites and correlation with experiments. J. Phys. Chem. B 115 (2011) 13320-13327.

59.

S. Parimal, S.M. Cramer, S. Garde, Application of a Spherical Harmonics Expansion Approach for Calculating Ligand Density Distributions around Proteins. J. Phys. Chem. B (2014).

60.

J. Woo, S. Parimal, M.R. Brown, R. Heden, S.M. Cramer, The effect of geometrical presentation of multimodal cation-exchange ligands on selective recognition of hydrophobic regions on protein surfaces. J. Chromatogr. A 1412 (2015) 33-42.

61.

J.A. Woo, H. Chen, M.A. Snyder, Y. Chai, R.G. Frost, S.M. Cramer, Defining the property space for chromatographic ligands from a homologous series of mixedmode ligands. J. Chromatogr. A 1407 (2015) 58-68.

62.

N. Chennamsetty, V. Voynov, V. Kayser, B. Helk, B.L. Trout, Prediction of aggregation prone regions of therapeutic proteins. J. Phys. Chem. B 114 (2010) 6614-6624.

63.

K. Srinivasan, S. Parimal, M.M. Lopez, S.A. McCallum, S.M. Cramer, Investigation into the Molecular and Thermodynamic Basis of Protein Interactions in Multimodal Chromatography Using Functionalized Nanoparticles. Langmuir 30 (2014) 13205-13216.

64.

H.S. Karkov, B.O. Krogh, J. Woo, S. Parimal, H. Ahmadian, S.M. Cramer, Investigation of protein selectivity in multimodal chromatography using in silico designed Fab fragment variants. Biotechnol. Bioeng. 112 (2015) 2305-2315.

65.

A. Hirano, T. Arakawa, T. Kameda, Interaction of arginine with Capto MMC in multimodal chromatography. J. Chromatogr. A 1338 (2014) 58-66.

66.

F. Dismer, J. Hubbuch, 3D structure-based protein retention prediction for ionexchange chromatography. J. Chromatogr. A 1217 (2010) 1343-1353.

30

67.

F.K. Insaidoo, M.A. Rauscher, S.J. Smithline, N.C. Kaarsholm, B.P. Feuston, A.D. Ortigosa, T.O. Linden, D.J. Roush, Targeted purification development enabled by computational biophysical modeling. Biotechnol. Prog. 31 (2015) 154-164.

68.

F. Fogolari, A. Brigo, H. Molinari, The Poisson–Boltzmann equation for biomolecular electrostatics: a tool for structural biology. J. Mol. Recognit. 15 (2002) 377-392.

69.

R. Kumari, R. Kumar, A. Lynn, g_mmpbsa- A GROMACS Tool for HighThroughput MM-PBSA Calculations. J. Chem. Inf. Model. 54 (2014) 1951-1962.

70.

P.A. Kollman, I. Massova, C. Reyes, B. Kuhn, S. Huo, L. Chong, M. Lee, T. Lee, Y. Duan, W. Wang, Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 33 (2000) 889-897.

71.

J.N. Israelachvili (2011) Intermolecular and surface forces (Academic press).

72.

T. Hou, J. Wang, Y. Li, W. Wang, Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations. J. Chem. Inf. Model. 51 (2010) 69-82.

73.

T. Hou, J. Wang, Y. Li, W. Wang, Assessing the performance of the molecular mechanics/Poisson Boltzmann surface area and molecular mechanics/generalized Born surface area methods. II. The accuracy of ranking poses generated from docking. J. Comput. Chem. 32 (2011) 866-877.

74.

D.P. Oehme, R.T. Brownlee, D.J. Wilson, Effect of atomic charge, solvation, entropy, and ligand protonation state on MM‐PB (GB) SA binding energies of HIV protease. J. Comput. Chem. 33 (2012) 2566-2580.

75.

C.M. Roth, K.K. Unger, A.M. Lenhoff, Mechanistic model of retention in protein ion-exchange chromatography. J. Chromatogr. A 726 (1996) 45-56.

76.

N. Boardman, S. Partridge, Separation of neutral proteins on ion-exchange resins. Biochem. J. 59 (1955) 543.

77.

W. Kopaciewicz, M. Rounds, J. Fausnaugh, F. Regnier, Retention model for highperformance ion-exchange chromatography. J. Chromatogr. A 266 (1983) 3-21.

31

78.

E.S. Parente, D.B. Wetlaufer, Relationship between isocratic and gradient retention times in the high-performance ion-exchange chromatography of proteins: theory and experiment. J. Chromatogr. A 355 (1986) 29-40.

79.

J.C. Love, L.A. Estroff, J.K. Kriebel, R.G. Nuzzo, G.M. Whitesides, Selfassembled monolayers of thiolates on metals as a form of nanotechnology. Chem. Rev. 105 (2005) 1103-1170.

80.

M. Mondello, G.S. Grest, E.B. Webb III, P. Peczak, S.T. Milner, Dynamics of nalkanes: Comparison to Rouse model. arXiv preprint physics/9802022 (1998).

81.

W.D. Cornell, P. Cieplak, C.I. Bayly, I.R. Gould, K.M. Merz, D.M. Ferguson, D.C. Spellmeyer, T. Fox, J.W. Caldwell, P.A. Kollman, A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 117 (1995) 5179-5197.

82.

J. Woo (2014) Engineering Biomolecular Interactions to Enhance Selectivity in Chromatographic Separations. (RENSSELAER POLYTECHNIC INSTITUTE).

83.

A.J. Patel, P. Varilly, S.N. Jamadagni, H. Acharya, S. Garde, D. Chandler, Extended surfaces modulate hydrophobic interactions of neighboring solutes. Proc. Natl. Acad. Sci. 108 (2011) 17678-17683.

84.

G.B. Sigal, M. Mrksich, G.M. Whitesides, Effect of surface wettability on the adsorption of proteins and detergents. J. Am. Chem. Soc. 120 (1998) 3464-3473.

85.

M. Rostkowski, M.H. Olsson, C.R. Søndergaard, J.H. Jensen, Graphical analysis of pH-dependent properties of proteins predicted using PROPKA. BMC Struct. Biol. 11 (2011) 1.

86.

H.J. Berendsen, D. van der Spoel, R. van Drunen, GROMACS: A message-passing parallel molecular dynamics implementation. Comput. Phys. Commun. 91 (1995) 43-56.

87.

B. Hess, C. Kutzner, D. Van Der Spoel, E. Lindahl, GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 4 (2008) 435-447.

88.

S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M.R. Shirts, J.C. Smith, P.M. Kasson, D. van der Spoel, GROMACS 4.5: a high-throughput and

32

highly parallel open source molecular simulation toolkit. Bioinformatics (2013) 29 (2013) 845-854. 89.

D.J. Evans, B.L. Holian, The nose–hoover thermostat. J. Chem. Phys. 83 (1985) 4069-4074.

90.

S. Nosé, A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 81 (1984) 511-519.

91.

M. Parrinello, A. Rahman, Crystal structure and pair potentials: A moleculardynamics study. Phys. Rev. Lett. 45 (1980) 1196.

92.

U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, A smooth particle mesh Ewald method. J. Chem. Phys. 103 (1995) 8577-8593.

93.

G. Hummer, S. Garde, A.E. Garcia, A. Pohorille, L.R. Pratt, An information theory model of hydrophobic interactions. Proc. Natl. Acad. Sci. U.S.A. 93 (1996) 89518955.

33

LIST OF FIGURES Figure 1. (A) Capto MMC. Interacting groups: phenyl group (circled green), amide linkage (circled blue), negatively charged carboxylic acid group (circled red) and thioalkyl tail (circled orange). (B) Nuvia cPrime. Interacting groups: phenyl group (circled green), amide linkage (circled blue), negatively charged carboxylic acid group (circled red). Cyan box shows the point of immobilization of these ligands to the base matrix. Figure 2. Flowchart showing method to generate bound structures of proteins to MM surfaces. Figure 3. Protein-Capto MMC SAM binding free energy G  ,   maps for four proteins. pH 5.0, NaCl concentration = Clow = 0.25M.  ,   are represented in units of radians. Figure 4. Linear relationship between calculated overall free energy and logarithm of experimental isocratic retention factors at 4 salt concentrations for (A) Capto MMC (R2=0.91) and (B) Nuvia cPrime (R2=0.93). pH 5.0. Free energies for each protein are evaluated at NaCl concentrations of 0.6M, 0.65M, 0.7M and 0.75M for Capto MMC and at 0.3M, 0.35M, 0.4M and 0.45M for Nuvia cPrime. Figure 5. Evolution of binding free energies of the “strongly bound” structures (filled circles) with increasing NaCl salt concentration. These “strongly bound” structures are the orientations that exhibited favorable free energy of binding at Clow of 0.25M. Proteins shown are Chymotrypsinogen-A and Horse Cytochrome C. Surface: Capto MMC, pH 5.0.  ,   are represented in units of radians. Figure 6. Predicted and experimental gradient elution salt concentrations for Ubiquitin and Ribonuclease A on (A) Capto MMC, (B) Nuvia cPrime. pH 5.0, salt: NaCl. Figure 7. Binding of Chymotrypsinogen-A to MM surfaces at pH 5.0. (A) Binding patches of various regions on protein surface to Capto MMC surface, (B) Binding patches of various regions on protein surface to Nuvia cPrime surface, (C) Electrostatic potential map, (D) Spatial aggregation propensity (hydrophobicity) map.

34

Figure 8. Binding of protein Horse Cytochrome C to MM surfaces at pH 5.0. (A) Binding patches of various regions on protein surface to Capto MMC surface, (B) Binding patches of various regions on protein surface to Nuvia cPrime surface, (C) Electrostatic potential map, (D) Spatial aggregation propensity (hydrophobicity) map. Figure 9. Binding of protein Ribonuclease-A to MM surfaces at pH 5.0. (A) Binding patches of various regions on protein surface to Capto MMC surface, (B) Binding patches of various regions on protein surface to Nuvia cPrime surface, (C) Electrostatic potential map, (D) Spatial aggregation propensity (hydrophobicity) map. Figure 10. Binding of protein Ubiquitin to MM surfaces at pH 5.0. (A) Binding patches of various regions on protein surface to Capto MMC surface, (B) Binding patches of various regions on protein surface to Nuvia cPrime surface, (C) Electrostatic potential map, (D) Spatial aggregation propensity (hydrophobicity) map.

35

Figr-1

Figure 1. (A) Capto MMC. Interacting groups: phenyl group (circled green), amide linkage (circled blue), negatively charged carboxylic acid group (circled red) and thioalkyl tail (circled orange). (B) Nuvia cPrime. Interacting groups: phenyl group (circled green), amide linkage (circled blue), negatively charged carboxylic acid group (circled red). Cyan box shows the point of immobilization of these ligands to the base matrix.

36

Figr-2

Figure 10. Binding of protein Ubiquitin to MM surfaces at pH 5.0. (A) Binding patches of various regions on protein surface to Capto MMC surface, (B) Binding patches of various regions on protein surface to Nuvia cPrime surface, (C) Electrostatic potential map, (D) Spatial aggregation propensity (hydrophobicity) map.

37

Figr-3

Figure 2. Flowchart showing method to generate bound structures of proteins to MM surfaces.

38

Figr-4

Figure 3. Protein-Capto MMC SAM binding free energy G  ,   maps for four proteins. pH 5.0, NaCl concentration = Clow = 0.25M.  ,   are represented in units of radians.

39

Figr-5

Figure 4. Linear relationship between calculated overall free energy and logarithm of experimental isocratic retention factors at 4 salt concentrations for (A) Capto MMC (R2=0.91) and (B) Nuvia cPrime (R2=0.93). pH 5.0. Free energies for each protein are evaluated at NaCl concentrations of 0.6M, 0.65M, 0.7M and 0.75M for Capto MMC and at 0.3M, 0.35M, 0.4M and 0.45M for Nuvia cPrime.

40

Figr-6

Figure 5. Evolution of binding free energies of the “strongly bound” structures (filled circles) with increasing NaCl salt concentration. These “strongly bound” structures are the orientations that exhibited favorable free energy of binding at Clow of 0.25M. Proteins shown are Chymotrypsinogen-A and Horse Cytochrome C. Surface: Capto MMC, pH 5.0.  ,   are represented in units of radians.

41

Figr-7

Figure 6. Predicted and experimental gradient elution salt concentrations for Ubiquitin and Ribonuclease A on (A) Capto MMC, (B) Nuvia cPrime. pH 5.0, salt: NaCl.

42

Figr-8

Figure 7. Binding of Chymotrypsinogen-A to MM surfaces at pH 5.0. (A) Binding patches of various regions on protein surface to Capto MMC surface, (B) Binding patches of various regions on protein surface to Nuvia cPrime surface, (C) Electrostatic potential map, (D) Spatial aggregation propensity (hydrophobicity) map.

43

Figr-9

Figure 8. Binding of protein Horse Cytochrome C to MM surfaces at pH 5.0. (A) Binding patches of various regions on protein surface to Capto MMC surface, (B) Binding patches of various regions on protein surface to Nuvia cPrime surface, (C) Electrostatic potential map, (D) Spatial aggregation propensity (hydrophobicity) map.

44

Figr-10

Figure 9. Binding of protein Ribonuclease-A to MM surfaces at pH 5.0. (A) Binding patches of various regions on protein surface to Capto MMC surface, (B) Binding patches of various regions on protein surface to Nuvia cPrime surface, (C) Electrostatic potential map, (D) Spatial aggregation propensity (hydrophobicity) map.

45