Molecular Simulations of the Pairwise Interaction of Monoclonal

Oct 28, 2014 - For a more comprehensive list of citations to this article, users are encouraged ... Structure and Relaxation in Solutions of Monoclona...
0 downloads 0 Views 7MB Size
Download PDF
Article pubs.acs.org/JPCB

Molecular Simulations of the Pairwise Interaction of Monoclonal Antibodies Mauro Lapelosa,† Thomas W. Patapoff,‡ and Isidro E. Zarraga*,† †

Department of Late Stage Pharmaceutical Development and ‡Department of Early Stage Pharmaceutical Development, Genentech Inc., member of Roche, South San Francisco, California 94080, United States S Supporting Information *

ABSTRACT: Molecular simulations are employed to compute the free energy of pairwise monoclonal antibodies (mAbs) association using a conformational sampling algorithm with a scoring function. The work reported here is aimed at investigating the mAb−mAb association driven by weak interactions with a computational method capable of predicting experimental observations of low binding affinity. The simulations are able to explore the free energy landscape. A steric interaction component, electrostatic interactions, and a nonpolar component of the free energy form the energy scoring function. Electrostatic interactions are calculated by solving the Poisson−Boltzmann equation. The nonpolar component is derived from the van der Waals interactions upon close contact of the protein surfaces. Two mAbs with similar IgG1 framework but with small sequence differences, mAb1 and mAb2, are considered for their different viscosity and propensity to form a weak interacting dimer. mAb1 presents favorable free energy of association at pH 6 with 15 mM of ion concentration reproducing experimental trends of high viscosity and dimer formation at high concentration. Free energy landscape and minimum free energy configurations of the dimer, as well as the second virial coefficient (B22) values are calculated. The energy distributions for mAb1 are obtained, and the most probable configurations are seen to be consistent with experimental measurements. In contrast, mAb2 shows an unfavorable average free energy at the same buffer conditions due to poor electrostatic complementarity, and reversible dimer configurations with favorable free energy are found to be unlikely. Finally, the simulations of the mAb association dynamics provide insights on the self-association responsible for bulk solution behavior and aggregation, which are important to the processing and the quality of biopharmaceuticals.

■

concentration,3 in which multiple interactions can exist. From an experimental perspective, static light scattering (SLS) and dynamic light scattering (DLS) are able to probe weak interactions, which are difficult to determine by X-ray crystallography. A small-angle neutron scattering study was able to explain rheology behavior of monoclonal antibodies (mAbs) due to self-association and aggregates.4,5 Although NMR techniques were successful in measuring interactions of proteins and molecules of small size,6 mAbs due to their high molecular weight still represent a significant challenge. The role of long-range electrostatics has been considered critical to the function and structure of proteins,7 but it can also contribute significantly to the kinetics of the dimerization process.8 Protein net charges alone, however, do not correlate well with weak interaction parameters such as kD or B22;3 therefore, higher order electrostatic multipoles are important to consider as well. The correlation of charges on nearest neighbor atoms is insufficient to determine the overall charge complementarity in a series of protein−protein interfaces investigated.9 Instead, the

INTRODUCTION The desire to gain insight into the dimerization and aggregation process of proteins has driven substantial experimental studies and retains a significant interest in biomedical science and drug development. It is well-known that protein−protein association encompasses a wide range of interactions. Depending on the strength and type of interaction, one can form either a substantial amount of oligomeric states or allow monomers to predominantly remain. Liquid biopharmaceutical formulations commonly contain a small amount of dimers and aggregates (typically < 3%) even after the purification process. At higher protein concentration required for subcutaneous injections, the growth rate of these dimers and aggregates can be pronounced and thus can impact the quality of biopharmaceuticals.1,2 Protein association can involve strong interactions in the case of tightly bound protein complexes characterized by a high degree of geometric complementarity and by the presence of covalent, or hydrophobic interactions, or hydrogen-bonds contacts. However, there are also many situations where the protein self-association is governed by weak interactions. These can include electrostatic and hydrophobic contacts in the limited area of interaction. Although these interactions are weak they can have a significant impact on solution properties at high © XXXX American Chemical Society

Received: August 29, 2014 Revised: October 23, 2014

A

dx.doi.org/10.1021/jp508729z | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

the dynamics of a fully atomistic model. This is in contrast to coarse-grained techniques used to simulate the behavior of mAbs in solution,30 where close range interactions are approximated by mAb bead models to reduce the computational burden and provide an approximate and qualitative determination of solution behavior and trends. To mitigate the high computational burden, an efficient searching algorithm is used with calculations on a grid taking into account the orientational dependence of all the atoms pairwise interactions. Extensive conformational sampling is employed to determine the free energy landscape. The mAbs are investigated at pH 6 and different ionic concentrations. mAb1 at low ionic strength presents favorable association free energy in comparison with mAb2. The methodology is successful for the determination of the low minimum energy configurations of the dimer for mAb1.

electrostatic potentials at larger protein−protein interfaces of more than just a few surface amino acids are required. Consideration of both local charge complementarity and far field electrostatic interactions is important in determining the electrostatic potential at the binding interface. The mAbs (Mw ≈ 145 kD) investigated in this study belong to the IgG1 class which is the most prevalent immunoglobulin in serum. This IgG1 framework has been widely employed in therapeutic applications and is part of a rapidly growing class of pharmaceutical products.10,11 It is known that mAbs can form dimers under various stress conditions (e.g., thermal stress, low or high pH, process stresses) typically encountered in pharmaceutical process development.2,3,12−15 These dimers are formed by the interaction of fragment antigen-binding (Fab) arms in the mAbs or else Fab−Fc interactions. In this computational study, mAb1 and mAb2 are selected because they present very different physicochemical properties in terms of viscosity and propensity to form self-associated oligomers.16−21 Indeed neutron spin echo measurements reveal that mAb1 forms reversible clusters at high concentrations substantially more than mAb2.5 Since the Fc region of the two mAbs is identical, the difference in clustering behavior must be related to the differences in the complementarity determining region (CDR) and Fv region of the mAbs that create a large charge asymmetry for mAb1. Molecular modeling of weak interactions together with the determination of free energy landscape and dimer configurations represent a significant challenge even with powerful and parallel computational resources. Traditional computational methods based on docking,22,23 which succeed in the case of high steric complementarity at the interaction surface between the two proteins, perform poorly in the case of low affinity interactions. Therefore, simulations involving atomistic details beyond surface complementarity are needed for a more accurate determination of minimum free energy configurations for mAb−mAb interactions, including local electrostatic complementarity and the impact of surface flexibility. Conformational sampling in synergy with a free energy scoring function provides the basis to analyze the process of dimerization for the mAbs, and a more detailed understanding of the solution behavior of mAbs. The mAbs simulated here, mAb1 and mAb2, have been well studied experimentally and provide a good test for the validity of the computational methodology. To model intermolecular interactions, previous computational studies have used a hybrid potential as Hamaker/ Lennard-Jones,24 where the hybrid model comprises a longrange treatment used to describe interactions between atoms separated by more than 6 Å, and a short-range part used for all closer interactions. This same model has been later further parametrized in an attempt to reproduce the experiments on protein−protein interaction so that it was able to determine B22 values, which provide a good description of proteins in solution.25 Following these studies other calculations of weak protein−protein interactions have been pursued with more emphasis in atomic details, and ultimately also directed to the development of new energetic potentials.26−29 Our work actually follows these last developments including atomistic details in the model to treat correctly long-range electrostatic effects, which also make the calculations more demanding in terms of computational power. In particular, the high molecular weight of a mAb and its highly anisotropic shape compared to globular proteins like lysozyme, make it challenging to simulate

■

THEORY AND METHODS Energy Function. The prediction of the free energy of association is achieved using an energy function that includes free energy components for steric overlap, polar interactions, and nonpolar interactions. The energy function is based on the following equation: ΔGassoc = ΔGoverlap + ΔGnonpolar + ΔGpolar

(1)

The ΔGoverlap represents the steric overlap energy evaluated on all atoms, ΔGnonpolar is the energetic component generated by van der Waals contacts formed upon binding. This represents the Lennard-Jones potential. ΔGpolar represents the electrostatic energetic component. The approach used to determine ΔGpolar (eq 2) for the full atomistic solvated system was based on the Poisson− Boltzmann (PB) equation.31−33 ΔGpolar =

1 2

∑ qiΦi

(2)

where the summation is carried over all the partial charges qi, and the potential Φi is precalculated according to eq 3. ∇·[ε(r )∇ϕ(r )] − k 2(r ) sinh ϕ(r ) = −4πρ f (r )

(3)

where ε the dielectric constant, k the Debye−Hü ckel parameter, and ρf represents the charge density. All the terms in this equation are position-dependent. A quantity that can be calculated to measure the weak interacting forces is the second virial coefficient (B22). Although B22 can be measured for mixed protein solutions, it can be also calculated for single-component protein solutions describing the interaction between two identical macromolecules. For monoclonal antibodies, the orientationally averaged nature of B22 can reveal the anisotropic character of protein interactions that cannot be approximated by spherical models.5 Since interactions involve limited local areas on the protein surface, refined molecular models are needed to capture the interaction details. Computationally, B22 is calculated according to B22 =

−1 16π 2M w 2

∫ (exp(−ΔGassoc/kT ) − 1) dΩ

(4)

where k is the Boltzmann constant, T is the absolute temperature, and Mw is the protein molecular weight. dΩ indicates that the integral is over all possible orientations. The magnitude of B22 reflects the size and direction of deviations from ideality of the solution experimentally, with an ideal B

dx.doi.org/10.1021/jp508729z | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

solution obeying the van’t Hoff relation. A positive B22 generally indicates repulsive interactions while a negative B22 implies attractive interactions. Computational Protocol. The homology models of mAb1 and mAb2 were prepared using as a template the human IgG1 crystal structure (PDB ID: 1HZH).34 The starting models of Fab generated with MODELLER software35 were superimposed onto the Fab regions of 1HZH and the atoms were replaced to generate the full model of the mAb (Figure 1).

stage of equilibration and production dynamics once the system was solvated, the cutoff distance for calculating electrostatic interactions directly was set to 9 Å. Particle mesh Ewald summation was used to calculate electrostatic interactions beyond this distance. The cutoff for calculating van der Waals interactions was also set to 9 Å. The production run was 200 ns long. After the MD simulations the equilibrated structures were selected from the trajectory using a hierarchical clustering technique. RMS was applied as distance metric on the 20000 trajectory snapshots with the total number of clusters set to eight to maximize the separation ratio among the clusters. The representative conformations from the most populated clusters for mAb1 and mAb2, considered the most probable structures in solution, were selected for the PB calculations and the conformational sampling procedure. The adaptive Poisson−Boltzmann Solver version 1.331 was used for generating the electrostatic potential for the selected structures of mAb1 and mAb2. For APBS calculations, the PDB2PQR version 1.538 was used with the AMBER forcefield36 to generate the atomic charges and the PQR file from the PDB coordinates. The PROPKA program39 was used to predict the pKa values and the protonation states using a cutoff of 6 Å for residue deprotonation energy at pH 6.0 for the respective mAbs. PROPKA predicts empirically pKa values given the structural details for ionizable chemical groups using six adjustable parameters and taking into account the pKa shifts due to energy components related to desolvation, hydrogen bonding, and intraprotein energies. These effects are determined by the local conformation and chemical nature of residues surrounding the specific ionization site. The pHspecific PQR file was subsequently used to calculate the electrostatic surface charge distribution of the protein with a linearized form of the PB equation and cubic B-spline discretization of the charge distributions. PB calculations were performed at 300 K with a dielectric constant of 78.0 for water and 2.0 for the protein interior with a grid of 200 Å for each axis. The ion concentrations varied from 150 mM to 15 mM with ionic radius of 2.0 Å. The conformational search for the exploration of the 6dimensional space was carried out at pH 6 and two ionic concentration (15 mM and 150 mM) conditions based on the Euler−Rodriguez axis-angle convention. The Euler parameters are a,b,c, and d, the rotation angle ψ for the atom i of the protein sampling state k, and the rotation matrix A:

Figure 1. Structural models of the Fab of mAb1 (top) and mAb2 (bottom). Positive and negative residues are drawn in blue and red, respectively. mAb1 presents changes that create charge asymmetry in comparison with mAb2. VL is negative, and VH is neutral in mAb1.

Hydrogens were added using AMBER tools with the force field ff99SB36 followed by a minimization for 5000 steps using a conjugate gradient. The atomistic model for each mAb includes about 300.000 atoms with TIP3P model used for water molecules. A energy minimization on the entire system, each using 5000 iterations of steepest descent followed by 1000 iterations of conjugate gradient, were carried out with the restraints on the solute reduced in increments of 5 kcal/mol in each round from 25 kcal/mol in the first to 0 kcal/mol in the last. A 10 ns MD equilibration followed in which the temperature of the system was increased from 100 to 300 K. The density of the system increased from 0.8770 g/cm3 to 1.0002 g/cm3 during equilibration of the solvated system. Molecular dynamics with a full atomistic and explicitly solvated model using AMBER1137 was employed to determine the equilibrated structures. An initial MD equilibration of the solvent was carried out with all protein and carbohydrate atoms restrained in a 50 kcal/ mol energy. Also the system was heated from 100 to 300 K over a period of 1.5 ps, with pressure held at 1 atm, and taking 1 fs time steps for integration of atomic equations of motion. Production MD simulations with pressure and temperature held constant at 1 atm and 300 K were carried out from this point on. For this stage of simulation, the SHAKE algorithm was used, and the size of the time steps for integrating the atomic equations of motion was increased to 2 fs. For each

⎡ x ⎤k ⎡ x ⎤0 ⎢ y ⎥ = A (k )⎢ y ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ z ⎦i ⎣ z ⎦i

where ⎡ a2 + b2 − c 2 − d 2 2(bc + ad) ⎤ 2(bd − ac) ⎢ ⎥ ⎥ A = ⎢ 2(bc − ad) a 2 + c 2 − b2 − d 2 2(cd + ab) ⎢ ⎥ 2 2 2 2 ⎢⎣ 2(bd + ac) 2(cd − ab) a + d − b − c ⎥⎦

(5)

and a = cos θ/2, b = kx sin θ/2, c = ky sin θ/2, d = kz sin θ/2. k⃗ is a unit vector. The rigid-body conformational search on the grid with three rotational and three translation dimensions has been performed using parallel threads subdividing the large conformational space, and consequently reconstructed once the energy for each grid point was calculated. Principal component analysis (PCA)40 was applied to further decompose the 6-dimensional C

dx.doi.org/10.1021/jp508729z | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Significant differences arise in the CDR in which the secondary structure of amino acids changes in the region that is more exposed to the solvent. The representative structure of mAb1 coming from cluster one shows high structural similarities to the crystal structure of Fab (Figure 3). The C-α RMSD between the two sets of atoms in the Fab is 0.78 Å.

space to the most significant three principal components (PCs), which span about 60% of the total variance of the data. We chose a set of x, y, z, and Euler parameters, to span the 1.5 million conformations as points in the 6-dimensional structure space. The six coordinates system is a lower limit for the unique determination of the relative positions of the about 21 000 atoms of each mAb. The PCA method is able to further reduce these six coordinates to the three most significant ones. The protein representations were produced with PyMOL41 and VMD.42

■

RESULTS Molecular Dynamics. Eight molecular dynamics simulations for 200 ns were used to determine the most probable structure in solution for mAb1, and mAb2 starting from an initial structure obtained by homology modeling. The two mAbs show large conformational changes in solution mainly due the flexible hinge loop region connecting the Fabs and the Fc, and both mAbs significantly reorganize the conformation from the initial structural model. Hierarchical clustering has been used to define the representative structure of the mAb to use for the electrostatic calculations. The larger cluster for mAb1 includes about 20% of the total snapshots sampled in the simulations. Clusters 2, 3, 4, 5, 6, 7, and 8 yielded respectively 18, 17, 17, 15, 6, 5, and 3% of the total population. On the other side, the larger cluster in mAb2 includes about 50% of the snapshots. The other mAb2 clusters 2, 3, 4, 5, 6, 7, 8 have respectively 21, 7, 7, 7, 7, 1, and 0% of the total populations. The representative structures of the five largest clusters for mAb1 and the largest cluster for mAb2 as seen in Figure 2 are

Figure 3. Superimposition of the representative structure of cluster 1 for mAb1 (orange) and the Fab crystal structure of mAb1 (cyan) (PDB ID: 2XA8).

PB Calculations. mAbs are initially positioned on a 3D-grid to calculate the electrostatic component of the free energy by solving the Poisson−Boltzmann equation. It is noticeable that mAb1 and mAb2 show substantial differences in the electrostatic potential at different pH with an ion concentration of 15 mM (Figure 4). mAb1 at pH 5 presents a large positive region

Figure 4. Representation of the electrostatic potential field lines colored with blue (potential > 1 kT/e) and red (potential < −1 kT/e) for mAb1 (top) and mAb2 (bottom) at 15 mM ion concentration. The panel on the left has pH 5, the middle one pH 6, and the right one has pH 7.

of the potential (blue) and a small negative region (red). Instead, the map changes at pH 6 where the negative region is more pronounced with a distinct separation. At pH 7, the positive region is less extended. For mAb2 at pH 5 the map presents mainly a large positive region. At pH 6, the Fc presents some red areas with extended intensity that increase at pH 7. Overall the electrostatics of mAb1 and mAb2 diverges significantly. Conformational Sampling. The conformational sampling simulations are performed for mAb1, and mAb2 at pH 6, chosen to compare the data with previous experiments used for these mAbs.17,19 Once the precalculated electrostatic potential is mapped on the 3D-grid, the conformational search is carried

Figure 2. Representative structures of mAb1 (top) and mAb2 (bottom) obtained from hierarchical clustering of MD trajectories.

used for the electrostatic calculations and the conformational sampling. Both the representative structures for mAb1 and mAb2 are overall more compact than the initial model, and the Fab arms extend themselves into the solvent perpendicularly to the Fc. The representative structure of mAb1 is overall similar to mAb2, although some differences are noticeable in the conformation and angles of Fabs relative to the Fc region. D

dx.doi.org/10.1021/jp508729z | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

out considering all the six degrees of freedoms (see Theory and Methods). The computational method developed allows the evaluation of the electrostatic potential without resolving the PB equation for each protein−protein conformation which is rather unfeasible computationally. The electrostatic interaction energy of a second protein is calculated by summing the contributions from each of its partial charges, in which each charge makes a contribution of partial charge on the second protein according eq 2, and is the electrostatic potential at that charge generated by the first protein. All the partial charges are considered to provide high resolution rather than just using the charges on ionizable residues as done previously.26 The methodological setup allows the generation of 1.5 million dimer conformations in each run. This specific approach saves a considerable amount of time and is possible only because the conformational flexibility is not allowed in either protein. Approximately thirty percent out of the total number of conformations for the mAbs are overlapping configurations and contribute to the excluded volume effect. By analogy van der Waals energies are evaluated on the 3D-grid around each protein to determine the nonpolar contribution to the association free energy, ΔGnonpolar. The behavior at pH 6 was simulated and compared with known rheological measurements of these mAbs used in pharmaceutical developments.16,17,19,21 The calculated values of ΔG̅ assoc obtained with pH 6 are in Table 1. mAb1 presents a

Figure 5. Binding energies distribution for mAb1 and mAb2 at 15 mM and 150 mM ion concentration.

Table 1. Average Association Free Energy for mAb1, and mAb2 at pH 6 Using Different Ion Concentrationsa ΔG̅ assoc

15 mMb

150 mMb

mAb1 mAb2

−0.75 ± 0.020 0.13 ± 0.014

−0.06 ± 0.018 0.02 ± 0.012

a

The error is the standard error of the mean. The average association free energy for mAb1 is calculated including representative structures from the five largest clusters. bValues in kcal/mol.

value of −0.75 kcal/mol as average interaction free energy at pH 6 with 15 mM of salt concentration. Changing the ion concentration to 150 mM brings the value of the association free energy down to −0.06 kcal/mol. mAb2 shows a positive value of 0.13 kcal/mol at 15 mM, which decreases to 0.02 kcal/ mol at 150 mM. Overall, from a quantitative perspective mAb1 presents more propensity to associate than mAb2 with a clear distinction between high and low ion concentrations. The two mAbs present a small difference at high ion concentration. These results show consistency with the electrostatic potential map shown previously in Figure 4. A free energy landscape can be reconstructed from the conformational search algorithm. Distributions of binding energies show the values of mAb1 and mAb2 to be close to 0 at high ion concentration (Figure 5). Instead, mAb1 binding energy distribution at 15 mM is shifted toward favorable free energy (lower than 0), and mAb2 shows unfavorable binding free energy. Four lowest minimum free energy dimer configurations for mAb1 are found (Figure 6). Fab−Fab interactions contribute for the most part to the favorable association free energy. The four configurations show a free energy of association of −0.9 kcal/mol (1), −1.4 kcal/mol (2), −1.2 kcal/mol (3), and −2.6 kcal/mol (4). In these configurations (1−2), the negative region of the electrostatic potential interacts with the positive region in a unique geometrical arrangement that can involve alternatively both

Figure 6. Four lowest free energy minima for mAb1 dimer association at ph6 with 15 mM ion concentration. Fab−Fab dimer configuration (top panel). Alternative Fab−Fab configuration (second middle panel). Configuration in which Fabs of the two mAbs are parallel to each other (third middle panel). Configuration double Fc−Fab mode of interaction (bottom panel).

Fab arms. This creates a dimer configuration with a radius of gyration of about 81 Å and 82 Å for the two Fab−Fab configurations of the dimer in Figure 6, a smaller value of about 70 Å coming from the measurement of a small-angle neutron E

dx.doi.org/10.1021/jp508729z | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

scattering characterization of the same mAbs.43 An alternative configuration of the mAb1 dimer presents the two molecules laying parallel on the z-axis displaced such that the positive and negative region of the Fab arm can interact with each other with a radius of gyration of about 55 Å. This configuration minimizes the surface exposure of the Fc arm, rich in hydrophobic residues, and the structure is more compact. The fourth configuration of the dimer presents a double Fc− Fab association with a radius of gyration of 71 Å. Figure 7 shows the field lines of the electrostatic potential for a dimer configuration of the mAb1 calculated solving the PB calculation

at pH 6 on the minimum free energy configuration of the dimer obtained with the conformational search algorithm. It is noticeable that the positive region and the negative one of the potential interacts for the two Fab arms. On the other side, the Fc region of each mAb stays with a strong positive region of the potential, and the small negative areas are present in the Fab on the opposite side from the binding region. Overall the dimer as a unit presents a large volume dominated by positive electrostatic potential with the exception of the two small negative regions. The field lines of the electrostatic potential calculated at pH 7 show a clear difference from the ones at pH 6. The extended negative area in the Fab disfavors the encounter of the two molecules in this conformational arrangement. The B22 calculations show a binding ability for mAb1 with a negative value of −9.26 × 10−05 at 15 mM ion concentration and −1.59 × 10−05 at 150 mM. These results are in overall agreement with previous B22 measurements obtained with dynamic light scattering techniques.18,21 The minimum free energy pathway for mAb1 is found starting from the unbound conformation using steepest descent (Figure 8). This shows how the structure reorients at medium range to assume the bound state.

Figure 8. Representative trajectory of minimum free energy pathway using the steepest decent approach from the unbound configuration to the bound state of the dimer of mAb1, extrapolated from the free energy landscape. One mAb molecule (receptor) involved in the conformational search is drawn in green and the binder is in cyan.

PCA method for mAb1 at 15 mM presents PCA1, which is mostly correlated with the z-axis of rotation and PCA2, which is mostly correlated with the x-rotation (Figure 9). PCA3 is mostly correlated with the energy. The histogram of PCA3 values shows a bimodal distribution behavior with most of the values around −1.5 and 1.5. The values of PCA1 around the zaxis are associated with the Fab−Fab orientation (Table 2).

■

DISCUSSION The overall results presented here show the feasibility of calculating minimum free energy configurations in the case of weak association between molecules, which is typically more difficult to perform than strongly bound proteins. The two mAbs studied here have a similar conformational arrangement, with minor differences in the relative position of Fab and Fc arms, and the most probable conformations were used in the simulation of pairwise interactions. The effects of the different CDR and Fv sequences were shown to have a substantial impact on the association properties of mAb1 and mAb2. While equilibrium conformations in solution help to provide a more accurate determination of electrostatic potential in comparison to those based on

Figure 7. Fields line representation of the electrostatic potential for one minimum free energy configuration of the dimer of mAb1 at pH 6 (top) and at pH 7 (bottom) and 15 mM of ion concentration. The blue color represents the positive region and the red one is the negative region of the potential. The threshold between −1 kT/e and +1 kT/e separates the positive and negative areas of the potential. F

dx.doi.org/10.1021/jp508729z | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

to pH 6 as seen in Figure 4. The change of protonation state in the pH shift modulates significantly the electrostatic potential in the CDR reducing the negative area of the potential in the Fab, so that the electrostatic complementarity increases. The conformational search results show the propensity to form a reversible dimer for mAb1 at low ion concentration as seen in previous experiments.17,19−21,45 The association is driven by the deep negative region of the electrostatic potential in the Fab region of mAb1 generated by the negative residues, and appears to result from a local dipole effect between the Fab arms. A Fab−Fc type of interaction is also obtained from the simulations and is related to charge−charge interaction. mAb1 has a theoretical charge of +23 at pH 6 with 15 mM of ion concentration, but it still interacts to form a dimer. mAb2 shows no favorable free energy for both ion concentrations tested. mAb2 shows a very positive theoretical charge at pH 6 with an ion concentration of 15 mM (+32). The overall theoretical charge does not seem to explain the dimerization and viscosity behavior for these mAbs. Instead, the separation and the spatial distribution of the electrostatic potential shed insight on the dimerization behavior. The self-association property of mAb1 in comparison with mAb2 relies on the spatial distribution of the electrostatic potential and the possibility of geometrical fit in the interaction. mAb2 does not self-interact mostly because it is difficult to find on average conformations, in which the potential is complementary on the interacting surfaces of the mAb. The change in ionic strength also affects the degree of interaction, in particular in the case of mAb1 the association becomes very weak at 150 mM demonstrating how high ionic strength can easily screen electrostatic interactions. This can produce interactions that are stabilizing, highly directional, and distant-dependent. It has been known that the polar nature of protein−protein interfaces allows greater electrostatic stabilization of complexes than is possible, for example, in the interior of proteins. This stabilization can result in part from the surface charge distribution exhibiting complementary electrostatic multipoles on opposing surfaces, though sufficiently weak to allow reversible association of the interacting macromolecules. The Fab−Fab interaction represents the principal mode of mAb1 self-association. This confirms the experimental findings where viscosity increases at high concentration were not reduced simply by the removal of the Fc region.45 The interactions observed between positive and negative residues along the Fab confirm previous evidence that suggested that the dipole−dipole interaction dominates the driving forces for the viscosity behavior of mAb1.46 In addition, the Fab−Fc mode of interaction also contributes to the favorable part of the free energy as it can be noticed by one of the minimum free energy configuration found in the simulations, and this is consistent with previous simulations performed with a low resolution coarse grain model.30 Also, Fc−Fc interactions are mostly unfavorable because a positive potential characterizes largely the Fc region in which a large positive patch is present. The explanation of prevalent Fab−Fab interaction for mAb1 relies on the evidence that when the two binding partners interact in this orientation they have more possibilities geometrically to interact than Fab−Fc. One of the minimum energy configurations of the dimer shows a radius of gyration of 71 Å, which is a relatively extended dimer configuration consistent with the experimental value of 69 Å at 10 mg/mL obtained using neutron spin echo.43

Figure 9. PCA analysis plots of mAb1 at 15 mM. PCA1 histogram is on the top panel, and PCA2 histogram is on the middle panel. PCA3 histogram is represented on the bottom panel.

Table 2. Canonical Correlation Analysis on PCAs variables

PCA1

PCA2

PCA3

1 2 3 4 5 6 7 8

−0.0847 −0.0026 −0.0039 0.0188 −0.3188 −0.4887 0.8075 −0.0028

−0.0000 −0.0000 −0.0000 0.0000 0.7512 −0.6527 −0.0985 −0.0000

0.6987 0.0932 0.1431 −0.6833 −0.0235 −0.0360 0.0595 0.1016

homology modeling due to solvent effects and flexibility,44 the slight conformational differences between the two mAbs considered here appear to play a minor role in the association properties of mAb1 and mAb2. In contrast, for mAb1 the histidines in the CDR-H3 play a crucial role in shifting the positive area of the potential in CDR loop 3 going from pH 7 G

dx.doi.org/10.1021/jp508729z | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Yearly et al.5 analyzed the protein−protein association behavior of mAb1 as function of the concentration using small angle neutron scattering (SANS).5 Overall, in their work it was found that the effective interaction of mAb1 at low concentrations was attractive, but this was confounded by repulsive interactions at high concentration. Fixing the theoretical charge number at +16, the net attraction strength decreases for mAb1 as the concentration increases. Therefore, for mAb1 proteins, the effective attraction strongly depends on its concentration no matter which charge number is used. If the experimental charge number can be considered to be a good approximation of the real charge number, the estimated net attraction strength based on the experimental charge number becomes very small at 150 mg/mL (the volume fraction is about 0.216), indicating that there is only a weak attraction between mAbs at this concentration. Our results can explain these trends, indeed the extended positive potential displayed by the dimeric unit mAb1 structure differs from the overall electrostatic field generated by the mAb1 monomer. At low concentration when mAb1 is mostly a monomer, and the potential is more heterogeneous (positive and negative regions). Once the mAbs organize themselves as a dimer unit the attraction becomes less pronounced as shown in Figure 7. The change of the potential from the monomer to the dimer unit may explain the reason of the decrease of attraction going from low to high concentration of mAb1. When the reversible dimers form, particularly for the ones based on Fab−Fab interactions with extended geometrical arrangements, the internal dipoles within the dimer nearly cancel out, leaving a strong monopole−monopole repulsion between dimers. This is evident in the experimental structure factor S(Q) obtained via SANS by Yearley et al. (2013).5 At 50 mg/mL, the mAb1 compressibility S(0) is substantially larger than the hard sphere model indicative of attraction, but for 100 and 150 mg/mL the S(Q) behavior suggests repulsion as the equilibrium shifts. Furthermore, the S(0) for mAb1 is about twice that of mAb2 for the 100 and 150 mg/mL cases, suggesting the freely diffusing mAb1 colloids (dimers) are about two times larger than mAb2 colloids (monomers), which was subsequently confirmed by neutron spin echo.43 Moreover, the electrostatic potential of the dimer unit suggests possible modes of interaction to form transient oligomers when mAbs are crowded at high concentration, indeed the presence of a small fraction of oligomers has been suggested by a SLS study.18

possible orientational configurations for the dimer. It was observed that the extended dimer configuration had complementary internal dipoles that roughly canceled, and the configuration is consistent with SANS measurements indicating dimer−dimer repulsion at high concentrations. The calculations have been successful to understand the reason behind high viscosity and propensity of self-association. Although other methodologies with different energy functions can be applied to reach the same goal.25,26,47−49 The analysis of binding energy distributions successful in protein−ligand binding problems50,51 provided very useful insights in the binding process for the weak protein−protein interaction such as this case. Minimum free energy pathways can be constructed and provide a picture of the orientation of the molecule through the binding process connecting the intermediate states and determining the kinetic obstacles to reach the final bound state. The PCA analysis suggested the importance of the Fab−Fab orientation, and was able to reproduce a previous study that used 2D vibrational correlation spectroscopy and PCA to show that the conformational heterogeneity was influenced by the concentration effect for mAb1.52 The B22 calculations were carried out successfully and provided an interesting comparison with experiments in solution of mAbs. The results followed previous B22 theoretical calculations with proteins of medium size,25−27,53 but to our knowledge this represents the first attempt to determine B22 for proteins of a large size such as mAbs. The B22 provided a distinction between mAb1 and mAb2. Also, B22 values bridge nicely with experiments in which it is not possible to measure the equilibrium between monomer and dimer or oligomers at proper experimental concentrations. Our results reproduce previous experimental measurements for mAb1,18,21 where B22 was obtained from static light scattering data. The B22 value for mAb1 was about 10−5 mol mL/g2 which is close to the calculated value in Table 3. The B22 value for mAb2 of 15 × Table 3. Calculated B22 Values That Provide a Measure of Protein−Protein Interactions 15 mMa

B22

−9.26 × 10 2.03 × 10−05

mAb1 mAb2 a

150 mMa −05

−1.59 × 10−05 3.15 × 10−06

Values in mol mL/g2.

10−5 mol mL/g2 is more positive than our value. Overall, the calculations showed predictive power for B22, and differences are expected because of the complexity of the calculations. Our model presents limitations similarly to other methods in terms of probing the internal flexibility of the proteins during the conformational search even though a long MD was run to determine the equilibrated structure. Allowing flexibility would create prohibitive calculations, which are already very expensive, but it can be critical with other proteins that undergo side chain rearrangements upon binding. The methodology developed showed results consistent with rheology, size-exclusion chromatography, and SANS experiments, and it is useful in providing insights in the design and development of therapeutical proteins. Furthermore, the computational study provides a good step forward in understanding detailed molecular mechanisms and statistical thermodynamics associated with the dimerization process. Mitigating dimer and aggregate formation is a critical aspect in designing high concentration formulations. Hence, the ability

■

CONCLUSIONS Two mAbs were computationally investigated in this work, and one mAb was found to have a higher propensity to form reversible dimers based on weak interactions, particularly local electrostatic dipole interactions involving the Fab arm. Our method successfully analyzed the molecular interaction in the association process using a procedure that couples free energy scoring function with conformational search. The computational methodology demonstrated how the association free energy differentiates mAbs in the dimerization process. Sampling conformations is shown to be key in binding problems. mAb1 and mAb2 show very different long-range electrostatic properties that drive the formation of the reversible dimer for mAb1, but not for mAb2. mAb1 negative patches in the CDR are likely to interact with the positive regions explaining the reasons behind the dimer formation. The conformational search provides evidence for many different H

dx.doi.org/10.1021/jp508729z | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Recombinant Humanized Monoclonal Antibody to Vascular Endothelial Growth Factor. Biochemistry 1999, 38, 13960−13967. (14) Paul, R.; Graff-Meyer, A.; Stahlberg, H.; Lauer, M. E.; Rufer, A. C.; Beck, H.; Briguet, A.; Schnaible, V.; Buckel, T.; Boeckle, S. Structure and Function of Purified Monoclonal Antibody Dimers Induced by Different Stress Conditions. Pharm. Res. 2012, 29, 2047− 2059. (15) Zarraga, I. E.; Taing, R.; Zarzar, J.; Luoma, J.; Hsiung, J.; Patel, A.; Lim, F. J. High Shear Rheology and Anisotropy in Concentrated Solutions of Monoclonal Antibodies. J. Pharm. Sci. 2013, 102, 2538− 2549. (16) Lilyestrom, W. G.; Yadav, S.; Shire, S. J.; Scherer, T. M. Monoclonal Antibody Self-Association, Cluster Formation, and Rheology at High Concentrations. J. Phys. Chem. B 2013, 117, 6373−6384. (17) Liu, J.; Nguyen, M. D. H.; Andya, J. D.; Shire, S. J. Reversible Self-Association Increases the Viscosity of a Concentrated Monoclonal Antibody in Aqueous Solution. J. Pharm. Sci. 2005, 94, 1928−1940. (18) Scherer, T. M.; Liu, J.; Shire, S. J.; Minton, A. P. Intermolecular Interactions of IgG1 Monoclonal Antibodies at High Concentrations Characterized by Light Scattering. J. Phys. Chem. B 2010, 114, 12948− 12957. (19) Yadav, S.; Liu, J.; Shire, S. J.; Kalonia, D. S. Specific Interactions in High Concentration Antibody Solutions Resulting in High Viscosity. J. Pharm. Sci. 2010, 99, 1152−1168. (20) Yadav, S.; Shire, S. J.; Kalonia, D. S. Factors Affecting the Viscosity in High Concentration Solutions of Different Monoclonal Antibodies. J. Pharm. Sci. 2010, 99, 4812−4829. (21) Yadav, S.; Laue, T. M.; Kalonia, D. S.; Singh, S. N.; Shire, S. J. The Influence of Charge Distribution on Self-Association and Viscosity Behavior of Monoclonal Antibody Solutions. Mol. Pharm. 2012, 9, 791−802. (22) Mandell, J. G.; Roberts, V. A.; Pique, M. E.; Kotlovyi, V.; Mitchell, J. C.; Nelson, E.; Tsigelny, I.; Ten Eyck, L. F. Protein Docking Using Continuum Electrostatics and Geometric Fit. Protein Eng. 2001, 14, 105−113. (23) Smith, G. R.; Sternberg, M. J. Prediction of Protein−Protein Interactions by Docking Methods. Curr. Opin. Struct. Biol. 2002, 12, 28−35. (24) Neal, B. L.; Asthagiri, D.; Lenhoff, A. M. Molecular Origins of Osmotic Second Virial Coefficients of Proteins. Biophys. J. 1998, 75, 2469−2477. (25) Asthagiri, D.; Neal, B. L.; Lenhoff, A. M. Calculation of ShortRange Interactions between Proteins. Biophys. Chem. 1999, 78, 219− 231. (26) Elcock, A. H.; McCammon, J. A. Calculation of Weak Protein− Protein Interactions: The pH Dependence of the Second Virial Coefficient. Biophys. J. 2001, 80, 613−625. (27) Elcock, A. H.; Sept, D.; McCammon, J. A. Computer Simulation of Protein−Protein Interactions. J. Phys. Chem. B 2001, 105, 1504− 1518. (28) Madura, J. D.; Briggs, J. M.; Wade, R. C.; Davis, M. E.; Luty, B. A.; Ilin, A.; Antosiewicz, J.; Gilson, M. K.; Bagheri, B.; Scott, L. R.; McCammon, J. A. Electrostatics and Diffusion of Molecules in Solution: Simulations with the University of Houston Brownian Dynamics Program. Comput. Phys. Commun. 1995, 91, 57−95. (29) Stark, A. C.; Andrews, C. T.; Elcock, A. H. Toward Optimized Potential Functions for Protein−Protein Interactions in Aqueous Solutions: Osmotic Second Virial Coefficient Calculations Using the MARTINI Coarse-Grained Force Field. J. Chem. Theory Comput. 2013, 9, 4176−4185. (30) Chaudhri, A.; Zarraga, I. E.; Kamerzell, T. J.; Brandt, J. P.; Patapoff, T. W.; Shire, S. J.; Voth, G. A. Coarse-Grained Modeling of the Self-Association of Therapeutic Monoclonal Antibodies. J. Phys. Chem. B 2012, 116, 8045−8057. (31) Baker, N. A. Poisson−Boltzmann Methods for Biomolecular Electrostatics. Methods Enzymol. 2004, 383, 94−118.

to predict computationally the binding free energy and the association behavior of mAbs in solution can be a powerful tool in helping design mAbs appropriate for high concentration formulations.

■

ASSOCIATED CONTENT

* Supporting Information S

Further computational details for the dipole moments of the minimum energy dimer conformations for mAb1, and further information on the clustering analysis are provided. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We would like to acknowledge Jacob Corn and Barthélemy Demeule for useful suggestions, and Yun Liu from the University of Delaware for discussions on the relevant neutron scattering data.

■

REFERENCES

(1) Harris, R. J.; Shire, S. J.; Winter, C. Commercial Manufacturing Scale Formulation and Analytical Characterization of Therapeutic Recombinant Antibodies. Drug Dev. Res. 2004, 61, 137−154. (2) Shire, S. J.; Shahrokh, Z.; Liu, J. Challenges in the Development of High Protein Concentration Formulations. J. Pharm. Sci. 2004, 93, 1390−1402. (3) Connolly, B. D.; Petry, C.; Yadav, S.; Demeule, B.; Ciaccio, N.; Moore, J. M. R.; Shire, S. J.; Gokarn, Y. R. Weak Interactions Govern the Viscosity of Concentrated Antibody Solutions: High-Throughput Analysis Using the Diffusion Interaction Parameter. Biophys. J. 2012, 103, 69−78. (4) Castellanos, M.; Pathak, J.; Leach, W.; Bishop, S.; Colby, R. Explaining the Non-Newtonian Character of Aggregating Monoclonal Antibody Solutions Using Small-Angle Neutron Scattering. Biophys. J. 2014, 107, 469−476. (5) Yearley, E. J.; Zarraga, I. E.; Shire, S. J.; Scherer, T. M.; Gokarn, Y.; Wagner, N. J.; Liu, Y. Small-Angle Neutron Scattering Characterization of Monoclonal Antibody Conformations and Interactions at High Concentrations. Biophys. J. 2013, 105, 720−731. (6) Jensen, M. R.; Ortega-Roldan, J.-L.; Salmon, L.; van Nuland, N.; Blackledge, M. Characterizing Weak Protein-Protein Complexes by NMR Residual Dipolar Couplings. Eur. Biophys. J. 2011, 40, 1371− 1381. (7) Gilson, M. K. Theory of Electrostatic Interactions in Macromolecules. Curr. Opin. Struct. Biol. 1995, 5, 216−223. (8) Janin, J. The Kinetics of Protein−Protein Recognition. Proteins 1997, 28, 153−161. (9) McCoy, A. J.; Chandana Epa, V.; Colman, P. M. Electrostatic Complementarity at Protein−Protein Interfaces. J. Mol. Biol. 1997, 268, 570−584. (10) Carter, P. J. Potent Antibody Therapeutics by Design. Nat. Rev. Immunol. 2006, 6, 343−357. (11) Carter, P. J.; Hazuda, D.; Wells, J. A. Next Generation Therapeutics. Curr. Opin. Chem. Biol. 2013, 17, 317−319. (12) Deperalta, G.; Alvarez, M.; Bechtel, C.; Dong, K.; McDonald, R.; Ling, V. Structural Analysis of a Therapeutic Monoclonal Antibody Dimer by Hydroxyl Radical Footprinting. MAbs 2013, 5, 86−101. (13) Moore, J. M.; Patapoff, T. W.; Cromwell, M. E. Kinetics and Thermodynamics of Dimer Formation and Dissociation for a I

dx.doi.org/10.1021/jp508729z | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

(32) Gilson, M. K.; Honig, B. Calculation of the Total Electrostatic Energy of a Macromolecular System: Solvation Energies, Binding Energies, and Conformational Analysis. Proteins 1988, 4, 7−18. (33) Nicholls, A.; Honig, B. A Rapid Finite Difference Algorithm, Utilizing Successive Over-Relaxation to Solve the Poisson−Boltzmann Equation. J. Comput. Chem. 1991, 12, 435−445. (34) Saphire, E. O.; Parren, P. W.; Pantophlet, R.; Zwick, M. B.; Morris, G. M.; Rudd, P. M.; Dwek, R. A.; Stanfield, R. L.; Burton, D. R.; Wilson, I. A. Crystal Structure of a Neutralizing Human IGG against HIV-1: a Template for Vaccine Design. Science 2001, 293, 1155−1159. (35) Eswar, N.; Webb, B.; Marti-Renom, M. A.; Madhusudhan, M. S.; Eramian, D.; Shen, M.-Y.; Pieper, U.; Sali, A. Comparative Protein Structure Modeling using Modeller. Curr. Protoc. Bioinf. 2006, Chapter 5, Unit 5.6. (36) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins 2006, 65, 712−725. (37) Case, D. A.; Darden, T. A.; Cheatham, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Crowley, M.; Walker, R. C.; Zhang, W.; et al., Amber 11. University of California: San Francisco, CA, 2010. (38) Dolinsky, T. J.; Czodrowski, P.; Li, H.; Nielsen, J. E.; Jensen, J. H.; Klebe, G.; Baker, N. A. PDB2PQR: Expanding and Upgrading Automated Preparation of Biomolecular Structures for Molecular Simulations. Nucleic Acids Res. 2007, 35, W522−W525. (39) Søndergaard, C. R.; Olsson, M. H. M.; Rostkowski, M.; Jensen, J. H. Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values. J. Chem. Theory Comput. 2011, 7, 2284−2295. (40) David, C. C.; Jacobs, D. J. Principal Component Analysis: a Method for Determining the Essential Dynamics of Proteins. Methods Mol. Biol. 2014, 1084, 193−226. (41) The PyMOL Molecular Graphics System, version 1.6; Schrödinger, LLC: New York, 2010. (42) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33−38. (43) Yearley, E.; Godfrin, P.; Perevozchikova, T.; Zhang, H.; Falus, P.; Porcar, L.; Nagao, M.; Curtis, J.; Gawande, P.; Taing, R.; Zarraga, I.; Wagner, N.; Liu, Y. Observation of Small Cluster Formation in Concentrated Monoclonal Antibody Solutions and Its Implications to Solution Viscosity. Biophys. J. 2014, 106, 1763−1770. (44) Chopra, G.; Summa, C. M.; Levitt, M. Solvent Dramatically Affects Protein Structure Refinement. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 20239−20244. (45) Kanai, S.; Liu, J.; Patapoff, T. W.; Shire, S. J. Reversible SelfAssociation of a Concentrated Monoclonal Antibody Solution Mediated by Fab−Fab Interaction That Impacts Solution Viscosity. J. Pharm. Sci. 2008, 97, 4219−4227. (46) Singh, S.; Yadav, S.; Shire, S.; Kalonia, D. Dipole−Dipole Interaction in Antibody Solutions: Correlation with Viscosity Behavior at High Concentration. Pharm. Res. 2014, 1−10. (47) Camacho, C. J.; Weng, Z.; Vajda, S.; DeLisi, C. Free Energy Landscapes of Encounter Complexes in Protein−Protein Association. Biophys. J. 1999, 76, 1166−1178. (48) Camacho, C. J.; Vajda, S. Protein Docking along Smooth Association Pathways. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 10636− 10641. (49) Kim, Y. C.; Tang, C.; Clore, G. M.; Hummer, G. Replica Exchange Simulations of Transient Encounter Complexes in Protein− Protein Association. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 12855− 12860. (50) Gallicchio, E.; Lapelosa, M.; Levy, R. M. The Binding Energy Distribution Analysis Method (BEDAM) for the Estimation of Protein-Ligand Binding Affinities. J. Chem. Theory Comput. 2010, 6, 2961−2977. (51) Lapelosa, M.; Gallicchio, E.; Levy, R. M. Conformational Transitions and Convergence of Absolute Binding Free Energy Calculations. J. Chem. Theory Comput. 2012, 8, 47−60.

(52) Kamerzell, T. J.; Kanai, S.; Liu, J.; Shire, S. J.; Wang, Y. J. Increasing IgG Concentration Modulates the Conformational Heterogeneity and Bonding Network That Influence Solution Properties. J. Phys. Chem. B 2009, 113, 6109−6118. (53) Velev, O. D.; Kaler, E. W.; Lenhoff, A. M. Protein Interactions in Solution Characterized by Light and Neutron Scattering: Comparison of Lysozyme and Chymotrypsinogen. Biophys. J. 1998, 75, 2682−2697.

J

dx.doi.org/10.1021/jp508729z | J. Phys. Chem. B XXXX, XXX, XXX−XXX