Quantifying the Risks of Asparagine Deamidation and Aspartate

Subscriber access provided by GAZI UNIV

Article

Quantifying Risks of Asparagine Deamidation and Aspartate Isomerization in Biopharmaceuticals by Computing Reaction Free Energy Surfaces Nikolay V. Plotnikov, Satish Kumar Singh, Jason C. Rouse, and Sandeep Kumar J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b11614 • Publication Date (Web): 04 Jan 2017 Downloaded from http://pubs.acs.org on January 5, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Quantifying Risks of Asparagine Deamidation and Aspartate Isomerization in Biopharmaceuticals by Computing Reaction Free Energy Surfaces Nikolay V. Plotnikov1,†, Satish Kumar Singh1, Jason C. Rouse2 and Sandeep Kumar1,* 1

Pharmaceutical Research and Development, Biotherapeutics Pharmaceutical Sciences, Pfizer Inc., 700 Chesterfield Pkwy West, Chesterfield, MO 63017, USA 2

Analytical Research and Development, Biotherapeutics Pharmaceutical Sciences, Pfizer Inc., 1 Burtt Rd., Andover, MA 01810, USA

†

Present address: R&D Platform Technology and Science, GlaxoSmithKline,1250 S Collegeville Rd, Collegeville, PA 19426 *Corresponding Author: [email protected] Tel: 636-329-2362

Plotnikov et al. Page 1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract

Early identification of asparagine deamidation and aspartate isomerization degradation sites can facilitate successful development of biopharmaceuticals. Several knowledge-based models have been proposed to assess these degradation risks. In this contribution, we propose a physics-based approach to identify the degradation sites based on the free energy barriers along the prechemical conformational step and along the chemical reaction path. These contributions are estimated from classical and QM/MM molecular dynamics simulations. Computed barriers are compared to those for reference reactions in water within GNG and GDG sequence motifs in peptides (which demonstrate the highest degradation rates). Two major factors decreasing the degradation rates relative to the reference reactions are: steric hindrance toward accessing reactive conformations and replacement of water by less polar side chains in the solvation shell of transition states. Among potential degradation sites in the complementarity determining region of trastuzumab and between two DK sites in glial cell derived neurotropic factor this method identified N30T, N55G, D102G, and D95K, respectively, in agreement with experiments. This approach can be incorporated in early computational screening of chemical degradation sites in biopharmaceuticals.


Page 2 of 40

Page 3 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Deamidation of asparagine (Asn) and isomerization of aspartate (Asp) are spontaneous chemical reactions involved in degradation of proteins. They are often observed during manufacturing, storage and shipping of biopharmaceuticals. These reactions can potentially decrease binding affinity and/or conformational stability1-6 of monoclonal antibodies (mAbs). They are also responsible for irreversible temperature-induced loss of activity of enzymes7. Typically, only some of Asn and Asp residues present in a protein are at high risk of degradation. As part of early-stage developability assessment, it is imperative to identify these high-risk degradation sites (or hotspots), especially when they are located at functional sites, such as complementarity determining regions (CDRs) of mAbs or at interfaces between protein domains. Ideally, this assessment would be performed at the in silico stage without the need to express the protein and experimentally assess the degradation pathway. Problematic Asn and Asp hotspots can be replaced by Gln and Glu8,9, respectively. However, these mutations might decrease binding affinity, and alternative approaches such as bulking up the n+1 residue10 have been explored. Experimental characterization of chemical degradation sites involves extensive forced degradation and peptide mapping studies,11 in which proteins are incubated under particular stress conditions (in a time-course format), digested enzymatically into small peptides and subsequently analyzed with liquid chromatography-mass-spectrometry. Samples of the progressively stressed protein are compared to the unstressed protein sample, which serves as the time-zero reference. Unfortunately, during the enzymatic digestion some sites undergo method-induced degradation12. Molecular modeling can help rationalize experimental observations and verify consistency of data. That is why computational tools are becoming integral in a comprehensive assessment of chemical stability of biopharmaceuticals13.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Oligopeptides in water have been extensively used as model systems for degradation of Asn and Asp residues in proteins.14,15 These studies have identified the size of the following (n+1) residue as the main risk factor for degradation. Peptides in which Asn is followed by glycine have the shortest half-lives on the order of a day. Histidine or serine in the n+1 position increases the half-life of Asn, τ1/2, up to ~10 days, while leucine prolongs it up to ∼100 days. These studies gave rise to sequence-based empirical rules for prediction of degradation hotspots in protein. However, such sequence-based rules significantly overestimate the risk of degradation. For example, 31 % of Asp hot spots and 43 % of Asn hot spots identified in mAbs with sequencebased approaches are false positives16. Another problem is “non-canonical” degradation hotspots in proteins (e.g. DK17, NN16) which deviate from sequence-based rules. The likely reason for these deviations is the protein tertiary structure, which can affect the degradation rates18,19 in a number of ways. To improve accuracy of risk assessment of chemical degradation, recently developed empirical models have been parametrized using a number of structure-based descriptors such as protein secondary structure, local flexibility and solvent exposure16,20,21. Some descriptors are based on factors affecting the degradation mechanisms. Interestingly, certain structural parameters that most closely describe chemical transformation are sometimes found to be non-relevant in the classification of hotspots14. While these empirical models are computationally affordable and often do not require detailed understanding of reaction mechanisms, the predictive power of these empirical approaches depends on the training set (its availability and quality) and on the ad hoc choice of model parameters. An alternative is to use physics-based approaches to estimate properties ab initio, when no immediately relevant stability data is available. Such approaches, however, are more


Page 4 of 40

Page 5 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


computationally expensive, require a detailed understanding (or exploration) of the reaction mechanism and availability of a high-resolution 3D structure. In return, these methods provide a great deal of theoretical insight into reaction mechanisms, which can be used to develop successful strategies for (subsequent) rational protein design and engineering. Several physics-based models of different complexity were proposed to study kinetics and mechanisms of chemical reactions involved in degradation22-25. However, most of these approaches did not explicitly take into account the effect of protein environment with a notable exception25 which still somewhat overlooked the effect of electrostatics on chemical reactions in proteins26,27. To the best of our knowledge, no physics-based algorithms are currently available to assess chemical degradation risks based on relevant free energy surfaces in proteins. To bridge this gap, we developed and explored applications of a model which estimates relevant free energy changes in proteins at an affordable computational cost. To achieve that, we utilize Molecular Mechanics (MM) and hybrid Quantum Mechanics/ Molecular Mechanics28 (QM/MM) free energy calculations to rank multiple degradation sites relative to GNG and GDG peptides in water, which demonstrate the fastest degradation rates. This approach identifies an experimentally known Asp isomerization hotspot between two potential degradation sites within DK sequence motifs in a therapeutic protein. It also correctly identifies Asn deamidation and Asp isomerization hotspots among potential degradation sites in various sequence motifs in the CDRs of trastuzumab. Finally, it reveals molecular origin of variable chemical stability of considered Asn and Asp residues. Methods Reaction Mechanism and the Rate Limiting Step



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Generally accepted reaction mechanism for deamidation of Asn and isomerization of Asp at physiological conditions is the succinimide-mediated path which proceeds via a nucleophilic attack29 (Figure 1), while other mechanisms, including direct hydrolysis of Asn and a nucleophilic attack of the backbone oxygen on the Cγ atom of Asn have been proposed16. At the first step of the succinimide-mediated mechanism, the proton of the backbone nitrogen is transferred to a general base, which makes the nitrogen a strong nucleophile. Once deprotonated, the nitrogen attacks the Cγ atom of Asn or Asp to yield a short-living tetrahedral intermediate. The tetrahedral intermediate then transforms into a stable succinimide intermediate, which eventually hydrolyzes into iso-Asp and/or Asp (Figure S1). The formation of the tetrahedral intermediate is believed to be the rate-limiting step at a near-neutral pH.23,30 In this work we consider a simplified succinimide-mediated mechanism25 which proceeds via the intramolecular proton transfer to the carbonyl oxygen of Asn or Asp (Figure 1). Based on intrinsic pKa values, the Asp sidechain can potentially act as a proton acceptor at slightly acidic pH. However, the most likely proton acceptors, particularly for Asn, are hydroxyl anions from the bulk solvent, even at low concentrations and at slightly acidic pH values (see SI for relevant estimates). The reaction rate of Asn deamidation depends on the solvent dielectric constant,31 which indicates that the reaction indeed involves heterolytic bond cleavage and a highly polar transition state. At the same time, the reaction rate of Asp isomerization shows weak or no dependence on the solvent dielectric constant31. In case of Asp isomerization, both final states and the TS are charged, and in case of Asn deamidation, the proton transfer reaction leads to charge separation. The Asn deamidation reaction is also accelerated at basic pH values, which indicates involvement of hydroxyl anions in the mechanism.


Page 6 of 40

Page 7 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Nucleophilic attack is most likely to occur when the sidechain adopts conformations that minimize the nucleophilic attack distance (dNuA) between the backbone nitrogen of the n+1 residue and the electrophilic Cγ atom (Figure 2). The dihedral angle χ comprised of C-Cα-CβCγ atoms is the most impactful on dNuA. This dihedral is related to the χ1 dihedral angle, conventionally defined as Nn-Cα- Cβ-Cγ. The dNuA distance also depends on the dihedral angle ψ comprised of N-Cα-C-Nn+1 atoms, where n is the number of Asn or Asp residue. In fact, the dihedral angle ψ′ comprised of Cβ-Cα-C-Nn+1 is more relevant to the orientation of the sidechain, however it is related to ψ via ~120˚ shift by molecular geometry. Thus, particular combinations of the two dihedral angles (χ and ψ) minimize the nucleophilic attack distance dNuA. The corresponding conformations are reactive conformations, since the chemical reaction is most likely to originate from them. Conformational flexibility, or absence of steric hindrance, determines ability of Asn (or Asp) to adopt a reactive conformation. The first reactive conformation at χ~60˚ (upper panels in Figure 2) corresponds to an extended backbone with ψ~180˚, while the second reactive conformation at χ~-60˚ corresponds to a compact backbone with ψ~0˚ (two bottom panels in Figure 2). The values of both dihedral angles are approximate estimates, which can be significantly shifted by protein environment. The protein scaffold dictates the relative free energies of conformations as well as the corresponding conformational barriers. Therefore, the energetics of accessing reactive conformations in a protein has to be taken in consideration when estimating the energetics for the rate-limiting step of degradation reactions. Conformational barriers or high free energies of reactive conformation (∆GCONF) can slow down the chemical reactions involved in degradation. Protein environment also interacts with the changing electronic density in the course of the Plotnikov et al. Page 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 40

reaction. The solvent dipoles have to be re-oriented in order to solvate the new charge distribution in the transition state(s). Re-orientation of water dipoles costs significantly less energy than in protein, where orientation of dipoles is restricted by the protein fold. This results in the reaction rates being faster in aqueous environments than in proteins (Figure 3). Degradation Reaction Rates Estimated from Transition State Theory Transition state theory32 predicts the rate constant of a chemical reaction from the activation free energies of the rate-limiting step. The expression for the rate constant can be written as: ݇=ߢ

݇஻ ܶ Δ݃‡ ሺߦሻ ∙ ݁‫ ݌ݔ‬ቆ− ቇ ℎ ݇஻ ܶ

(1)

Where k is the reaction rate constant and ∆g‡ is its activation energy, which is a difference between the free energy level of the transition state and of the reactants state along a process coordinate ξ. The transmission coefficient κ is a constant which depends on the reaction type and on choice of the reaction coordinate; kB is the Boltzmann constant; T is the thermodynamic temperature and h is the Planck constant. Eq. 1 can be rearranged to express the ratio of two rate constants (or half-lives) for a reaction in different environments, e.g. in water and in protein. This ratio, to a good approximation, is given by the change in the activation free energy, ∆∆g‡: ‡ ݇௉ோை் Δ݃‡ − Δ݃ௐ஺் ΔΔ݃‡ = ݁‫ ݌ݔ‬ቆ− ௉ோை் ቇ = ݁‫ ݌ݔ‬ቆ− ቇ ݇ௐ஺் ݇஻ ܶ ݇஻ ܶ

(2)

Reactions of Asn deamidation and Asp isomerization proceed most rapidly within GNG and GDG motifs in small peptides in water. Therefore, we can take experimental degradation rates and calculated activation energies in water as reference points (that is, kWAT is the rate constant in a peptide in water). The smaller the difference in activation energies between a reaction in protein and the reference reaction in water, the faster the degradation reaction proceeds.


Page 9 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Furthermore, the calculated change in the activation free energy can provide an estimate for the ratio of half-lives (see SI for details). Free Energies of Reactive Conformations Computed from Potential of Mean Force Some Asn and Asp residues in a protein crystal structure are already found in reactive conformations, while others need to cross a conformational barrier to access a reactive conformation. Reactive conformations can also have a higher free energy, than non-reactive ones. Potential degradation sites are initially examined with MOE201533 by performing a rigid potential energy scan along the χ dihedral angle of Asn and Asp sidechain. Next, we compute probability densities for ψ and χ dihedral angles (Figures S2-S4) from a molecular dynamics trajectory in VMD34 or cpptraj.35 This analysis of accessibility of reactive conformations is limited by presence of high conformational barriers in proteins, which make sampling with regular molecular dynamics inefficient. To improve the conformational sampling, we employ the umbrella sampling36 approach. That is, a harmonic bias is added to the original molecular mechanics potential, EMM: ‫ܧ = ܧ‬ெெ + ‫ܭ‬ଵ ሺψ − ψ଴ ሻଶ + ‫ܭ‬ଶ ሺχ − χ଴ ሻଶ

(3)

We choose χ0 and ψ0 close to the values observed in the starting structure in order to avoid high energy gradients. Then we increase χ0 from the original value (in the relaxed structure) by increment of 5˚ and -5˚ to cover the range from 0 to 360 degrees. We use K1=200 kcal/(mol·rad) and K2=0 since rotation around the ψ angle is generally associated with higher energy barriers as it significantly disturbs the protein fold. The resulting molecular dynamics trajectories are propagated sequentially with GPU-accelerated37 Amber 1438. Once completed, the potential of mean force along the χ is constructed using the WHAM39 free energy estimator (while also monitoring changes in ψ). Plotnikov et al. Page 9 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 40

Reaction Free Energy Surfaces from QM/MM Accessibility of reactive conformations does not automatically mean that a given Asn or Asp residue is at high degradation risk. It simply means that there is no steric hindrance for the residue to reach the reactants state of the chemical step. However, a high free energy penalty of reactive conformation (∆GCONF) means that the residue is stable to degradation via the studied succinimide-mediated pathway. Those residues which are in or can easily access reactive conformations (within a few thermal units of energy, kBT) are further examined to estimate the activation and the reaction free energies at the chemical step. These free energies are evaluated from the QM/MM

28,40

free energy profiles which are

computed for each potential degradation site. The QM/MM approach allows us to quantify the effect of solvation on changing electronic density of reacting species at an affordable computational cost. The free energy profiles are computed along the minimum free energy path which was identified on the 2D free energy surface (see Figure S5) with the nucleophilic attack dNuA and the proton transfers dPT distances as reaction coordinate25. This is accomplished by introducing a harmonic bias for both order parameters to the QM/MM potential. ‫ܧ = ܧ‬ொெ/ெெ + ‫ܭ‬ଷ ሺdே௨஺ − d଴ே௨஺ ሻଶ + ‫ܭ‬ସ ሺd௉் − d଴௉் ሻଶ

(4)

We compute the free energy along the vector comprised of dPT and dNuA which were obtained from the 2D free energy surface (Figure S5). This is done in order to reduce the computational cost. We use K3=K4=250 kcal/(mol·Å2) and a semi-empirical PM641 hamiltonian in Amber QM/MM implementation42. The free energy profiles are estimated with a series of free energy estimators using a protocol described elsewhere43. Additional computational details are given in the SI.


Page 11 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Oligopeptides in water containing Asn and Asp followed by glycine (GNG and GDG) exhibit the highest degradation rates15 with the lowest activation energies. That is why we choose these reactions as reference points for ranking chemical degradation sites in proteins. That is, we estimate the difference in the chemical barrier between potential degradation sites in protein and GNG (GDG) peptide in water. It is important, however, to use as similar computational protocol as possible (e.g. size of QM and MM regions) for all reactions, including the reference reaction in water, to minimize non-systematic errors and to maximize systematic error cancellation. Results Asp isomerization in human glial cell derived neurotropic factor At first, we tested this approach on a human glial cell derived neurotropic factor (GDNF) protein. It has a known Asp isomerization hotspot, Asp9517, which is found within an uncommon (based on sequence-based trends in oligopeptides15) D95K sequence motif. The GDNF protein also has another aspartyl residue, Asp80, within the D80K sequence motif (Figure 4A), which is chemically stable. To the best of our knowledge, a crystal structure is not available for a human GDNF homodimer. We generated a homology-based structural model for the human GDNF using a hybrid template built by superposing a human monomer (PDB ID 2V5E44) on a rat homodimer (PDB ID 1AGQ45), which has 94 % pairwise identity. The crystal structure for the human monomer has a partially resolved Lys sidechain at the D95K site, and only one out of four monomers of the rat GDNF found in the crystal asymmetric unit has a resolved D95K site. Despite a very high degree of identity between the rat and the human GDNF, the initial conformations of the loops containing D95K sites are quite different, while D80K sites which lie within the helices are essentially structurally identical (Table 1). The D95K site in both homologs



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 40

is originally in a non-reactive conformation. However the rigid potential energy scan along the χ dihedral angle of Asp showed that a reactive conformation corresponds to a minimum of energy. The D80K site was observed only in a non-reactive conformation. Thus, the conformational flexibility of Asp and Lys sidechains is very different between D80K and D95K sites. D95K sites have several minima on the rigid potential energy scan along the χ dihedral angle (corresponding to non-reactive and reactive conformations) while D80K sites have a single minimum (nonreactive conformation). The probability density maps for ψ and χ dihedral angles (obtained from molecular dynamics) also indicate that D80K sites have their conformational space limited to a single non-reactive conformation (Figure S2) while D95K sites access both reactive and nonreactive conformations. Most importantly, the free energy of reactive conformation (∆GCONF) is ~1 kcal/mol for D95K and ~ 10 kcal/mol for D80K (Figure 5A). An order of magnitude difference in the free energy results in a dramatic difference in chemical stability between these two sites. Thus, we estimated the chemical barrier for the D95K site only and compared it to the barrier in Ace-GDG-Nme peptide in water. The free energy barrier in the peptide is ~5 kcal/mol lower compared to the D95K site in GDNF (Figure 6A and Table 3), which is in qualitative agreement with experimental estimates for half-lives (Table S1). Degradation sites in complementarity determining region of trastuzumab Next, we applied the proposed algorithm to characterize potential degradation sites for Asn and Asp in CDR loops of trastuzumab (which were defined according to the MOE internal annotation scheme). A crystal structure of the Fab region of trastuzumab in a complex with its antigen, HER2 receptor, is available (PDB ID: 1N8Z46). Degradation sites in trastuzumab have been experimentally characterized4,8,11,47. There are two Asn deamidation sites and one Asp


Page 13 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


isomerization site in the CDR loops: N30T motif in the CDR1 of the light chain (VL), N55G motif in the CDR2 of the heavy chain (VH) and D102G in the CDR3 of VH4 (Table S1). These sites are located at the binding interface between the Fab region and the HER2 receptor (Figure 4B). In addition to the experimentally determined degradation sites, other potential degradation sites in the CDR loops of trastuzumab include: D28V in CDR1 of VL, N28I and D31T in CDR1 of VH, D62S in CDR2 of VH and D108Y in CDR3 of VH. Structural parameters for these sites are given in Table 1. The primary degradation site, N30T, is in a reactive conformation, and a few chemically stable sites, such as D31T, are in non-reactive conformations. Curiously: a majority of chemically stable sites (N28I, D28V, D62S and D108Y) are also in reactive conformations with dNuA being significantly smaller or comparable to the original values (Table 1) in experimentally observed degradation sites such as D102G and N55G. Probability densities computed for ψ and χ dihedral angles from molecular dynamics (Figure S3-S4) also show that many of these sites remain in reactive conformations. This analysis also revealed differences in the conformational flexibility between these sites, with D102G being the most flexible one with several thermally accessible conformers. To determine the free energy of reactive conformation, ∆GCONF, we computed the free energy profiles along the χ dihedral angle (Figure 5B). The estimated free energies revealed that reactive conformations for N30T and for D102G are the most probable conformations, while for N55G and D62S reactive conformations are ~1 and ~2 kcal/mol higher. Finally, we computed the free energy profiles after bringing residues into corresponding reactive conformations (Figure 6B and Tables 2, S4). The resulting order of chemical barriers is consistent with experimental estimates for half-lives. The reaction free energy barrier (∆g‡) in D62S puts it among other degradation sites (D102G and D95K in GDNF), however, the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 40

conformational step, ∆GCONF, elevates the transition state by ~ 2 kcal/mol. Similarly, N55G and, surprisingly, N28I, have chemical barriers comparable to the N30T site. The reaction free energy change, ∆Gr, on the other hand, is the lowest for the N30T site, which indicates a strong driving force for the cyclization step. Thus, computed chemical barriers and reaction free energy changes distinguish experimentally known hotspots (Table S1) among potential degradation sites (Table S4). Despite the easy accessibility of reactive conformations, the chemical barriers for nondegradable residues are significantly higher compared to the reference reaction in water. In addition to chemical barriers (∆g‡) for experimentally known degradation sites being lower, the reaction free energies (∆Gr) for those sites are also lower (Tables 1 and S4) indicating a stronger driving force for the reaction. In other words, all hotspots show a greater stabilization of the transition state and of the product state. For the Asp isomerization reaction the endothermicity is consistently higher than for the Asn deamidation reaction, both for peptides in water and for proteins. The main difference between the formed tetrahedral intermediates is the net negative charge on the Asp-derived tetrahedral intermediate. The reaction free energy change for the Asp isomerization is highly affected by the difference in the solvation free energy between the reactants with the Asp side chain and the tetrahedral intermediate product. For Asn deamidation reactions, the free energy change is less positive, indicating a higher thermodynamic stability of the corresponding tetrahedral intermediate. Discussion Factors that Determine Reaction Rates: Steric Hindrance and Transition State Solvation Structure-based descriptors which correlate with high rates of chemical degradation16,20,21,48 can be classified into two categories: those describing conformational flexibility of degradation


Page 15 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


sites and those describing their solvent exposure. In our model, we distinguish a pre-chemical conformational step and a chemical step in the succinimide-mediated degradation mechanism. The steric hindrance toward adapting a reactive conformation affects the conformational step, and is captured via the free energy of reactive conformation, ∆GCONF, on the free energy surface along χ and ψ dihedral angles. Difference in solvation between the reactants and the transition state affects the chemical barrier, ∆g‡, while difference in solvation between the reactants and the products states affects the reaction free energy, ∆Gr. These parameters are estimated from the free energy profile along the reaction path. Steric hindrance49 restricts the conformational space of Asn and Asp residues (low-energy values of χ and ψ), and might prevent them from accessing reactive conformations (Figure 2). This can happen due to a high conformational barrier (kinetic control) or due to a high free energy difference between conformations (thermodynamic control). Steric hindrance is particularly manifested for residues with bulky sidechains and inside tightly packed protein cores, where rotation along some bonds can potentially lead to atomic clashes. Furthermore, charged sidechains might prefer certain rotamers within the protein framework50, which offer better solvation free energy gains due to preferential electrostatic interactions, such as saltbridges, or higher solvent-exposure. In other words, a chemical reaction might be sterically allowed or forbidden depending on the reactants ability to adopt a reactive conformation. We demonstrated that the D80K site in GDNF has a prohibitively high free energy penalty for accessing reactive conformations as compared to the D95K site. An order of magnitude difference in the conformational free energy change (∆GCONF) can be attributed to different secondary structures between D95K (loop) and D80K (helix). However, small but still significant free energy penalties are also observed for degradation sites in loops, for example, D62S and N55G sites in



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 40

trastuzumab. These free energy penalties elevate the transition state energy for the chemical step and, thus, decrease the reaction rates (Figure 3). Differences in solvation free energies between the reactant and the transition states affects reaction rates even if no steric hindrance prevents reactants from reaching reactive conformations. The difference in solvation is particularly manifested for reactions which proceed via heterolytic bond cleavage and involve charge separation. Asn deamidation is an example of such a reaction, where proton transfer leads to charge separation at the transition state. Enzymatic reactions in proteins proceed faster due to electrostatic stabilization of the transition state by evolutionarily pre-organized protein dipoles26. Degradation reactions, on the other hand, are non-enzymatic, and therefore the protein dipoles are not pre-organized to stabilize the transition state. Re-orienting protein dipoles to a new charge distribution at the transition state costs more energy than re-orienting water dipoles. Consequently, the activation barrier in protein is higher than in water. We chose the QM/MM truncation scheme which creates QM regions of the same size for various sequence motifs (Figures S6). We found that this is important after we examined the effect of the QM region size on absolute and relative values of QM/MM activation and reaction free energies (Figures S7, S8 and Table S4). The small QM region (CH3-CO-NHCH(CH2CONH2)-CO-NH-CH3) extends from the Cα atom of the n-1 residue up to the Cα atom of the n+1 residue (Figure S7-B). It includes Asn or Asp side chain, while all other side chains are described on the MM level. Thus, this QM/MM truncation scheme also produces QM regions independent on nature of flanking residues. One disadvantage of this scheme is that the methyl group (which results from the Cα atom of the n+1 residue) has a strong inductive effect, which increases pKa value of the backbone nitrogen. The small QM model gives exaggerated estimates


Page 17 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for the difference in the activation free energy relative to the reference reaction (Table S4). We found that extending the QM region up to the Cα atom of the n+2 residue (CH3-CO-NHCH(CH2CONH2)-CO-NH-CH2-CO-NH-CH3) significantly lowers the activation free energy and the reaction free energy. In earlier studies, the ψ2 dihedral angle (comprised of Nn+1-Cα-C-Nn+2 atoms) was suggested to contribute towards increased acidity of the backbone nitrogen22,25. Interestingly, our calculations showed no significant effect of ψ2 on the reaction free energy profiles for either reaction conformation (Table S4). Certain conformers of serine and threonine side chains were found to correlate with increased degradation risks18. While inclusion of side chains indeed lowers the barrier, it also makes the QM region dependent on the sequence motifs and complicates a systematic comparison across all sites. While this difficulty might be overcome by comparing barriers from different proteins for degradation sites within the same sequence motifs, we ensure the consistency by describing the n+1 sidechain on the MM level. In this way, the QM region for all potential degradation sites is comprised of the same atoms similar to the reference reactions within GNG or GDG peptides. Solvent accessibility and number of water molecules Given a random orientation of dipoles in water and in protein, water has a higher dielectric constant (ε=80) than protein (ε~4), even considering that the latter is ill-defined. Water dipoles should require less energy to be polarized to a new charge distribution at the transition state, than protein dipoles (Figure 3), which orientation is dictated by the protein structure. Therefore, solvent inaccessible Asn or Asp residues in proteins are less likely to undergo chemical degradation because this reaction proceeds through the charged transition states and intermediates (Figure S1), and the rate of Asn deamidation is dielectric dependent31. In addition, solvent-exposed residues are more likely to be conformationally flexible.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 40

Solvent exposure and number of water molecules in the 1st solvation shell of the reactive conformation are popular descriptors in predictive models for chemical degradation. Electrostatic interactions are long-ranged. Therefore, not only the 1st solvation shell (the amount of adjacent water molecules) affects the activation barrier, but also the total amount of water in the solvation shell of the transition state. The amount of water molecules in the solvation shell is the greatest for GNG (or GDG) peptides in water (Table 3). This can qualitatively explain why bulky side chains of the n+1 residue increase the activation energy compared to glycine. Furthermore, there is a strong correlation between the number of water molecules in the solvation sphere and the chemical barrier with the exception of D95K hotspot in GDNF which has a low number of water molecules in the 1st solvation shell. However, the positive charge on lysine side chain Nξ atom might, possibly, compensate for the lack of water molecules. The free energy barrier and the reaction free energy in GDNF are close to those for known Asp isomerization sites of trastuzumab (D102G) and are below those for stable sites in trastuzumab. In our model, QM/MM free energy surfaces quantitatively capture both the short- and the long-ranged electrostatic interactions. Algorithm quantifying degradation risks We examined the application of the proposed approach using two examples of biologics with determined crystal structures and previously characterized degradation sites. We were able to distinguish degradable and non-degradable sites in GDNF within the same DK sequence motif. Moreover, we were able to rank Asn and Asp residues in CDR loops of trastuzumab in agreement with experimental data. Due to its physics-based nature, this approach can be applied to assessment of chemical stability of Asn and of Asp in other amino acid-based modalities


Page 19 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


including therapeutic peptides, vaccines and engineered enzymes. A complete computational workflow is described below (Figure 7). We start from examining structural models of proteins. For potential hotspots we measure initial values of χ and ψ dihedral angles (and dNuA) to identify if a residue is in a reactive conformation (Figure 2). This step also includes a rigid potential energy scan along the χ dihedral angle for Asn or Asp sidechain. Initial values for potential degradation sites in GDNF and in trastuzumab are given in Table 1. It is hard, however, to quantify which residues are chemically unstable solely based on the initial assessment of accessibility of reactive conformations and on the distance for the nucleophilic attack, dNuA, even if this information is combined with the primarily structure. However, the analysis of the protein structure explains why it is problematic to predict chemical degradation risks using only sequence-based descriptors. During the initial assessment step, we also generate probability densities for χ and ψ dihedral angles from unbiased molecular dynamics with the explicit solvent to examine accessibility of reactive conformations (Figures S2-S4). In GDNF, only D95K sites were observed in reactive conformations (Figure S2), although non-reactive conformations had higher probability densities. In the CDRs of trastuzumab, almost all potential chemical degradation sites (except D62S) were found near reactive conformations (Figures S3 and S4). Here one can notice that for a few residues results obtained from the crystal structure with the rigid potential energy scan are different from the probability densities obtained from molecular dynamics. This can be attributed to, either or both, insufficient MD sampling, difference in solvation models and deficiency of energy minimization technique in proteins. Next, we compute the free energy cost of moving residues found in non-reactive conformations to reactive conformations. This is achieved by



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 40

constructing the free energy profile along the χ dihedral angle. For example, the free energy cost of accessing reactive conformation for D80K of GDNF is prohibitively large (Figure 5A), while for D62S and N55G of trastuzumab it is relatively moderate (Figure 5B). For residues which are either in reactive conformations or can access reactive conformations at a low free energy cost, we compute the free energy profiles along the chemical reaction path. Next, we compare the activation and the reaction free energies to values obtained for the reference reactions in GNG (or GDG) peptide in water, which exhibits the highest degradation rate (Table S1). This allows us to quantify the corresponding degradation risks. The free energy barrier height (∆g‡) along the chemical step is found to be the most important factor affecting the degradation rate (Table 3). The reaction free energy change (∆Gr) is the second important factor, which describes the driving force of a reaction (Table S4). For some residues, most notably for D62S, the free energy change at the conformational step (∆GCONF) seems to be critical for the residue chemical stability. Residues which are known to degrade experimentally have the smallest difference in both the activation free energy and the reaction free energy change as compared to the reference reaction for GNG (or GDG) peptide in water in addition to being in a reactive conformation. Limitations of the method Numerical accuracy of our approach is affected by several factors: limited knowledge and simplification of the reaction mechanism, empiricism in modeling MM and QM/MM potential energy surfaces, sampling and free energy estimation errors. To minimize the systematic error we apply an identical computational protocol for all degradation sites across different proteins and compare results to the reference reaction in water.


Page 21 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Quantitative description of chemical reactions is one of the key challenges in computational chemistry51. In this work we used a semi-empirical PM6 hamiltonian within the QM/MM scheme to describe the reaction potential energy surface. This level of theory allows for an extensive screening of multiple degradation sites at an acceptable computational cost. The difference between activation energies for the GNG and for the GDG peptide is significantly overestimated, since the NG peptide is known to degrade faster than the DG peptide. Based on difference in half-lives, the activation free energy difference for DG should be about 2 kcal/mol more than that for NG, which is around 20-22 kcal/mol31 (Table S1 and the SI). This quantitative discrepancy is most likely to be due to the reaction mechanism, particularly, in the approximation for the proton transfer path (compare Figure 1 and S1). Another limitation of the reaction mechanism depicted in Figure 1 is that the intramolecular proton transfer leads to a preferential configuration of a chiral center at Cγ, since rotation along Cβ - Cγ is frozen. The resulting diastereomer has a configuration, which depends on the starting reactive conformation. Thus, this algorithm can be refined by incorporating emerging mechanistic insight on the degradation reaction mechanism. The level of theory can be further improved either perturbatively43 or by performing electronic structure calculations on GPUs52. Convergence of the free energy calculations in proteins, particularly with explicit solvent, is another challenge51. Here, averaging over multiple replicas, extending simulation times and utilizing enhanced sampling algorithms can help minimize the associated errors. Some solutions employed in this work included: use of multiple starting configurations, GPU-accelerated molecular dynamics, and larger integration time step (via hydrogen mass repartitioning). Chemically degradable residues are often found in protein loops, which conformations can be difficult to predict correctly with the homology modeling and are often unresolved even in high-



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 40

resolution crystal structures53. This algorithm is not expected to be more sensitive to quality of homology models than previously proposed structure-based empirical algorithms, which is important since crystal structures for biologic drug candidates are rarely available at early discovery stages. The resulting error can be estimated by performing studies using multiple structural models. In our algorithm, a reactive conformation is the reactants state for the chemical step. The chemical step here is described by the minimum free energy path identified on the PM6/MM 2D free energy surface for Asp-isomerization at D95K modeled with the small QM region including the lysine sidechain. Our assumption is that both deamidation and degradation reactions at other sites in different proteins (described with QM regions of different sizes) follows this path closely. Some reactive conformations also have to be brought to the starting values at this path (dPT=1.36 Å and dNuA=3.11 Å) via energy minimization with the energy cost for this process neglected. This cost can be estimated via the restraint release approach or by computing the minimum energy path from every reactive conformation with the nudged elastic band or similar approach. Advantages of the algorithm The approach described in this contribution can help identify chemical degradation hotspots even if no experimental data is available. Despite a number of assumptions, our current model has identified Asn and Asp hotspots in semi-quantitative agreement with experimental data. It also revealed the molecular origin of main structural descriptors (size of flanking residues, flexibility of the degradation sites and their solvent accessibility) and provided valuable mechanistic insight on the succinimide-mediated degradation pathway. In order to form a cyclic intermediate, Asn and Asp sidechains should be able to adapt reactive conformations, which we


Page 23 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


describe by two dihedral angles. If rotation around these angles toward a reactive conformation, in particular around the χ, is limited by the steric hindrance, then the degradation pathway is disabled. On the other hand, certain conformations of the ψ angle (Figure 2) and its rigidity can also significantly decrease the degradation rate (e.g. N55G of trastuzumab). Only those residues which can reach reactive conformations proceed to the chemical step. Solvent accessibility affects the electrostatic stabilization of the transition state. The QM/MM activation and reaction free energies correlate with the number of water molecules in the solvation shell of reactive conformation for many degradation hotspots. When a part of the solvation shell is replaced by protein, the barrier increases. Bulky side chains of the n+1 residue displace water from the nearest solvation shell and increase the activation barriers. Conclusions In this contribution we presented a physics-based approach to predicting degradation sites in protein sequences which relies on modeling the actual processes of Asn deamidation and Asp isomerization directly in proteins. We demonstrated that conformational and chemical free energy barriers and free energy changes at the conformational and at the chemical steps can not only rationalize experimental data but also predict chemical degradation risks for Asn and Asp residues. Our findings are summarized in a computational algorithm, which, while still requires further improvements, can already be used to rank chemical liability of potential degradation sites of Asn deamidation and Asp deamidation in protein-based biologics. The method is general and applicable to any protein-based biopharmaceuticals, it doesn’t require any experimental information for parameterization and is ab initio in this sense. Calibration of the algorithm using existing experimental data can further improve accuracy of predictions. Supporting Information



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Additional details on developed computational protocols and on used experimental data are available as Supporting Information. This information is available free of charge via the Internet at http://pubs.acs.org Acknowledgement NP is grateful to the Pfizer postdoctoral program for financial support. We would like to thank our colleagues at Pfizer from Biotherapeutics Pharmaceutical Sciences and at Global Biotechnology Therapautics for critical discussions. Will Hickman, Nicholas Labello and Michael Miller and Pfizer HPC for their help with HPC resources and software.


Page 24 of 40

Page 25 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figures

Figure 1 Reaction scheme describing the succinimide-mediated reaction mechanism of Asndeamidation and of Asp-isomerization. Asn residue (n) is shown in blue color, while the n+1 residue is shown in red color, and n-1 and n+2 residues are shown in black. The proton of the backbone nitrogen is transferred to the carbonyl oxygen of the amide group of Asn (or to the carboxylate group of Asp). The proton transfer is followed by the nucleophilic attack of the deprotonated nitrogen on the Cγ atom of the Asn/Asp sidechain. Newly forming chemical bonds are shown with red dash lines.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 40

Figure 2 Newman projections describing two idealized reactive conformations and showing the dihedral angles that minimize the distance for the nucleophilic attack (dNuA) of the backbone nitrogen (in red) on Cγ of Asn or Asp sidechain. The first reactive conformation (upper two structures) corresponds to the sidechain dihedral angle χ~60˚ ( it is comprised of C-Cα-CβCγ atoms and is shown in the upper left Newman projection) while the backbone dihedral angle ψ~180˚ (this angle is comprised of Nn-Cα-C- Nn+1 atoms, shown in the upper right projection). The second reactive conformation (the two lower structures) adopts χ ~ -60˚ (the lower left projection) at ψ ~ 0˚ (the lower right projection).


Page 27 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3 Effects of solvent and of steric hindrance on the reaction rate of chemical degradation. (A) Polar water solvates an ionic transition state better than a less polar protein environment, which results in a higher free energy barrier in protein (dashed blue line) compared to the reference reaction in water (red line). Additionally, water dipoles are easier to re-orient to the new charge distribution in the transition state compared to protein dipoles, which orientation is determined by the protein structural context. (B) Steric hindrance at the pre-chemical conformational step increases the free energy of reactive conformation. Greater steric hindrance (dashed blue line) results in the overall free energy barriers being elevated higher compared to a residue with lower steric hindrance (red line) even though the activation free energies at the chemical step are the same.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 40

Figure 4 (A) Human glial cell-derived neurotropic factor is a homodimer, which has an Asp isomerization site D95K (red bold font, underlined in the sequence). The same molecule has another D80K motif within its sequence, which is chemically stable (shown in black bold font, underlined in the sequence). (B) Hypervariable loops of complementarity determining region (CDR) of trastuzumab contain a number of potential degradation sites, among which only three are experimentally known to degrade (Asn deamidation in N30T of the light chain, N55G of the heavy chain, and Asp isomerization in D102G of the heavy chain). These sites (shown with balls and sticks) are at the interface with the antigen, HER2 receptor (depicted with electrostatic surface, where “+” is blue and “-“ is red). Amino acid sequences for the variable region of the light chain (VL) and of the heavy chain (VH) are given below the figure, with CDR loops of VL shown in purple and CDR loops of VH shown in tan, except the H3 loop which is shown in red.


Page 29 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5 Free energy profiles quantify steric hindrance which is encountered by Asn and Asp sidechains as they rotate along Cα-Cβ bond (χ dihedral angle) while moving into reactive conformations. Shaded areas of the same color indicate error bars. (A) Average over monomers A and B of human GDNF homodimer for D95K hot spot (red circles) and for D80K site (green circles). The reactive conformation for D95K is around 60˚ (ψ~180˚) and is significantly lower than the one for stable D80K sites (~300˚ for ψ~0˚). (B) Trastuzumab: N30T of VL (red solid triangles) and D102G of VH (green circles) exist in reactive conformations at χ~60˚ and ψ~180˚. N55G of VH (black flipped triangles) requires ~0.8 kcal/mol to access the reactive conformation at χ~-90˚ and ψ~10˚, while D62S of VH (blue circles) requires ~ 1.8 kcal/mol to reach the reactive conformation at χ~-70˚ and ψ~-20 ˚.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 40

Figure 6 QM/MM free energy profiles along the chemical reaction path from corresponding reactive conformations (zero levels) for Asn deamidation (A) and for Asp isomerization (B). (A) Asn deamidation: in GNG peptide in water (solid black line and squares, average over two reactive conformations for ψ2~0° and for ψ2~180°); in the CDR of trastuzumab: N55G site of VH (solid blue line with squares), N30T of VL (solid red line with squares) and N28I of VH (solid green line with squares). (B) Asp isomerization: in GDG peptide in water (solid black with empty circles); in the CDR of trastuzumab: D102G VH (solid wine red line with empty triangles), D62S VH (solid orange line with empty flipped triangles), D31T of VH (solid green line), D28V of VL (dashed green line) and D108Y of VH (solid blue line with empty rhombs); in GDNF: D95K (solid red line with empty squares).


Page 31 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7 The proposed algorithm for prediction of chemical degradation sites using the physicsbased approach. This algorithm includes two major steps to assess chemical degradation risks of Asn deamidation and Asp isomerization. First, using available structural model we determine which residues are in reactive conformations from which chemical reaction can proceed. For the residues in non-reactive conformations we compute the free energy of moving them to the reactive conformation, ∆GCONF, taking into account protein environment. If a residue is either in a reactive conformation or can access one at low free energy cost comparable to the reference value in water (∆GCONF_WAT), we compute the QM/MM free energy profile along the reaction path and compare the activation (∆g‡) and the reaction free energies (∆Gr) in protein to values for the reference reaction in GNG (GDG) peptide in water (∆g‡WAT and ∆Gr_WAT), which exhibits the highest degradation rate. If those values in protein are comparable to the reference values in Plotnikov et al. Page 31 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

water, a residue is at high degradation risk. If values in protein are much higher, the residue is at low risk.


Page 32 of 40

Page 33 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Tables Table 1 Values of conformational and chemical order parameters for trastuzumab from PDB ID 1N8Z and for human GDNF from PDB ID 2V5E. Values reported after the back slash indicate the minimum energy value obtained using rigid potential energy scan along the χ dihedral angle of the Asn/Asp sidechain (with Amber10 force field as implemented in MOE 2015). Sequence motif ݀ே௨஺ / Å ߯° ߰° ߰ᇱ° 28 N I CDR1VH -175/71-74 132 -102 4.734/3.479 31 D T CDR1VH -168 1 129 4.912 N55G CDR2 VH -87 8 136 4.032 62 D S CDR2 VH 149/-68 -21 101 4.724/3.304 D102G CDR3 VH 112/-76 -7 125 4.766/3.805 D108Y CDR3 VH -72 -26 99 3.405 D28V CDR1 VL 55 131 -101 3.079 N30T CDR1 VL 74 -114 9 3.091 80 D K GDNF 168 -36 88 4.690 D95K GDNF 159/53 147 -89 4.408/3.061 1 80 D K rat GDNF 168 -37 86 4.664 1 95 D K rat GDNF 163/68 -87 35 4.243/3.068 1 Using structure for rat GDNF PDB ID 1AGQ 94% sequence identity, 98 % similarity. 1.6 RMSD to human GDNF (PDB ID 1N8Z).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 40

Table 2 Summary of calculated data on activation (∆g‡) and reaction free energies for chemical (∆Gr) and for conformational step (∆GCONF). ∆g‡ ACE-GNG-NME 29.1 55 Trastuzumab N G.VH 35.1 Trastuzumab N30T.VL 36.4 Trastuzumab N28I.VL 37.7 Peptide GDG 21.4 Trastuzumab D102G.VH 27.3 95 GDNF D K 26.1 Trastuzumab D62S.VH 26.3 31 Trastuzumab D T.VH 33.0 28 Trastuzumab D V.VH 35.1 Trastuzumab D108Y.VH 33.8 1 values are given in kcal/mol

∆Gr 4.2 12.7 6.6 14.8 5.3 13.6 11.6 16.1 22.5 22.2 23.9


∆GCONF 0 0.8 0 0 0 0.9 1.8 -

Page 35 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3 Number of water molecules (average ± standard deviation) in the 1st and in the 2nd solvation shell of a reaction center in reactive conformations from QM/MM trajectories used to construct the free energy surface for the chemical step. ݊തௐ஺் < 3.4 Å1 Ace-GDG-Nme 7.0±1.4 1 102 D G CDR3 VH 6.9±1.4 D62S CDR2 VH 5.6±1.4 28 D V CDR1 VL 2.4±0.8 D31T CDR1 VH 3.9±1.3 108 D Y CDR3 VH 1.4±0.7 95 D K GDNF 3.3±1.2 Ace-GNG-Nme 4.2±1.5 30 N T CDR1 VL 3.3±1.3 N55G CDR2 VH 1.4±1.1 28 N I CDR1VH 1.1±1.0 1 Within the reaction center (Cγ and O of Asn or Asp and reactive conformation

݊തௐ஺் < 5.0 Å 23.6±2.4 23.5±2.4 19.1±2.9 8.6±1.6 11.4±2.3 5.5±1.7 12.6±2.4 23.2±2.5 16.1±2.6 13.0±2.3 9.5±2.0 N, H of the n+1 residue) in the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 40

References (1) Haberger, M.; Bomans, K.; Diepold, K.; Hook, M.; Gassner, J.; Schlothauer, T.; Zwick, A.; Spick, C.; Kepert, J. F.; Hienz, B.et al. Assessment of Chemical Modifications of Sites in the Cdrs of Recombinant Antibodies mAbs 2014, 6, 327-339. (2) Vlasak, J.; Bussat, M. C.; Wang, S.; Wagner-Rousset, E.; Schaefer, M.; KlinguerHamour, C.; Kirchmeier, M.; Corvaïa, N.; Ionescu, R.; Beck, A. Identification and Characterization of Asparagine Deamidation in the Light Chain Cdr1 of a Humanized Igg1 Antibody Analytical Biochemistry 2009, 392, 145-154. (3) Bults, P.; Bischoff, R.; Bakker, H.; Gietema, J. A.; van de Merbel, N. C. LcMs/Ms-Based Monitoring of in Vivo Protein Biotransformation: Quantitative Determination of Trastuzumab and Its Deamidation Products in Human Plasma Analytical Chemistry 2016, 88, 1871-1877. (4) Harris, R. J.; Kabakoff, B.; Macchi, F. D.; Shen, F. J.; Kwong, M.; Andya, J. D.; Shire, S. J.; Bjork, N.; Totpal, K.; Chen, A. B. Identification of Multiple Sources of Charge Heterogeneity in a Recombinant Antibody Journal of Chromatography B: Biomedical Sciences and Applications 2001, 752, 233-245. (5) Zhang, A.; Hu, P.; MacGregor, P.; Xue, Y.; Fan, H.; Suchecki, P.; Olszewski, L.; Liu, A. Understanding the Conformational Impact of Chemical Modifications on Monoclonal Antibodies with Diverse Sequence Variation Using Hydrogen/Deuterium Exchange Mass Spectrometry and Structural Modeling Analytical Chemistry 2014, 86, 3468-3475. (6) Yan, Y.; Wei, H.; Fu, Y.; Jusuf, S.; Zeng, M.; Ludwig, R.; Krystek, S. R.; Chen, G.; Tao, L.; Das, T. K. Isomerization and Oxidation in the Complementarity-Determining Regions of a Monoclonal Antibody: A Study of the Modification–Structure–Function Correlations by Hydrogen–Deuterium Exchange Mass Spectrometry Analytical Chemistry 2016, 88, 2041-2050. (7) Ahern, T. J.; Klibanov, A. M. The Mechanisms of Irreversible Enzyme Inactivation at 100c Science 1985, 228, 1280. (8) Yang, Y.; Zhao, J.; Geng, S.; Hou, C.; Li, X.; Lang, X.; Qiao, C.; Li, Y.; Feng, J.; Lv, M.et al. Improving Trastuzumab’s Stability Profile by Removing the Two Degradation Hotspots Journal of Pharmaceutical Sciences 2015, 104, 1960-1970. (9) Smith, R. C.; O'Bryan, L. M.; Mitchell, P. J.; Leung, D.; Ghanem, M.; Wilson, J. M.; Hanson, J. C.; Sossick, S.; Cooper, J.; Huang, L.et al. Increased Brain Bio-Distribution and Chemical Stability and Decreased Immunogenicity of an Engineered Variant of Gdnf Experimental Neurology 2015, 267, 165-176. (10) Patel, C. N.; Bauer, S. P.; Davies, J.; Durbin, J. D.; Shiyanova, T. L.; Zhang, K.; Tang, J. X. N+1 Engineering of an Aspartate Isomerization Hotspot in the ComplementarityDetermining Region of a Monoclonal Antibody Journal of Pharmaceutical Sciences 2016, 105, 512-518. (11) Chelius, D.; Rehder, D. S.; Bondarenko, P. V. Identification and Characterization of Deamidation Sites in the Conserved Regions of Human Immunoglobulin Gamma Antibodies Analytical Chemistry 2005, 77, 6004-6011. (12) Diepold, K.; Bomans, K.; Wiedmann, M.; Zimmermann, B.; Petzold, A.; Schlothauer, T.; Mueller, R.; Moritz, B.; Stracke, J. O.; Mølhøj, M.et al. Simultaneous Assessment of Asp Isomerization and Asn Deamidation in Recombinant Antibodies by Lc-Ms Following Incubation at Elevated Temperatures PLoS ONE 2012, 7, e30295. Plotnikov et al. Page 36 ACS Paragon Plus Environment

Page 37 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(13) Developability of Biotherapeutics Computational Approaches; Sandeep, K.; Singh, S., Eds.; CRC Press Taylor & Francis Group: Boca Raton, FL, 2016. (14) Robinson, A. B.; McKerrow, J. H.; Cary, P. Controlled Deamidation of Peptides and Proteins: An Experimental Hazard and a Possible Biological Timer Proceedings of the National Academy of Sciences 1970, 66, 753-757. (15) Robinson, N. E.; Robinson, A. B. Molecular Clocks Proceedings of the National Academy of Sciences 2001, 98, 944-949. (16) Sydow, J. F.; Lipsmeier, F.; Larraillet, V.; Hilger, M.; Mautz, B.; Mølhøj, M.; Kuentzer, J.; Klostermann, S.; Schoch, J.; Voelger, H. R.et al. Structure-Based Prediction of Asparagine and Aspartate Degradation Sites in Antibody Variable Regions PLoS ONE 2014, 9, e100736. (17) Hui, J. O.; Chow, D. T.; Markell, D.; Robinson, J. H.; Katta, V.; Nixon, L.; Chang, B. S.; Rohde, M. F.; Haniu, M. Identification of Asp95as the Site of Succinimide Formation in Recombinant Human Glial Cell Line-Derived Neurotrophic Factor Archives of Biochemistry and Biophysics 1998, 358, 377-384. (18) Kossiakoff, A. A. Tertiary Structure Is a Principal Determinant to Protein Deamidation Science 1988, 240, 191-194. (19) Stephenson, R. C.; Clarke, S. Succinimide Formation from Aspartyl and Asparaginyl Peptides as a Model for the Spontaneous Degradation of Proteins Journal of Biological Chemistry 1989, 264, 6164-6170. (20) Robinson, N. E. Protein Deamidation Proceedings of the National Academy of Sciences 2002, 99, 5283-5288. (21) Robinson, N. E.; Robinson, A. B. Deamidation of Human Proteins Proceedings of the National Academy of Sciences 2001, 98, 12409-12413. (22) Radkiewicz, J. L.; Zipse, H.; Clarke, S.; Houk, K. N. Neighboring Side Chain Effects on Asparaginyl and Aspartyl Degradation: An Ab Initio Study of the Relationship between Peptide Conformation and Backbone Nh Acidity Journal of the American Chemical Society 2001, 123, 3499-3506. (23) Peters, B.; Trout, B. L. Asparagine Deamidation: Ph-Dependent Mechanism from Density Functional Theory Biochemistry 2006, 45, 5384-5392. (24) Kaliman, I.; Nemukhin, A.; Varfolomeev, S. Free Energy Barriers for the NTerminal Asparagine to Succinimide Conversion: Quantum Molecular Dynamics Simulations for the Fully Solvated Model Journal of Chemical Theory and Computation 2010, 6, 184-189. (25) Ugur, I.; Marion, A.; Aviyente, V.; Monard, G. Why Does Asn71 Deamidate Faster Than Asn15 in the Enzyme Triosephosphate Isomerase? Answers from Microsecond Molecular Dynamics Simulation and Qm/Mm Free Energy Calculations Biochemistry 2015, 54, 1429-1439. (26) Warshel, A.; Aqvist, J.; Creighton, S. Enzymes Work by Solvation Substitution Rather Than by Desolvation Proc. Natl. Acad. Sci. U. S. A. 1989, 86, 5820-5824. (27) Warshel, A. Energetics of Enzyme Catalysis Proceedings of the National Academy of Sciences 1978, 75, 5250-5254. (28) Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme J. Mol. Biol. 1976, 103, 227-249.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 40

(29) Geiger, T.; Clarke, S. Deamidation, Isomerization, and Racemization at Asparaginyl and Aspartyl Residues in Peptides. Succinimide-Linked Reactions That Contribute to Protein Degradation Journal of Biological Chemistry 1987, 262, 785-794. (30) Capasso, S.; Mazzarella, L.; Sica, F.; Zagari, A.; Salvadori, S. Kinetics and Mechanism of Succinimide Ring Formation in the Deamidation Process of Asparagine Residues Journal of the Chemical Society, Perkin Transactions 2 1993, 679-682. (31) Brennan, T. V.; Clarke, S. Spontaneous Degradation of Polypeptides at Aspartyl and Asparaginyl Residues: Effects of the Solvent Dielectric Protein Science 1993, 2, 331-338. (32) Eyring, H. The Activated Complex and the Absolute Rate of Chemical Reactions Chemical Reviews 1935, 17, 65-77. (33) 2015 ed.; Chemical Computing Group Inc. : Montreal, QC, Canada, 2016. (34) Humphrey, W.; Dalke, A.; Schulten, K. Vmd: Visual Molecular Dynamics J. Mol. Graphics 1996, 14, 33-38. (35) Roe, D. R.; Cheatham, T. E. Ptraj and Cpptraj: Software for Processing and Analysis of Molecular Dynamics Trajectory Data Journal of Chemical Theory and Computation 2013, 9, 3084-3095. (36) Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling J. Comput. Phys. 1977, 23, 187-199. (37) Salomon-Ferrer, R.; Götz, A. W.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with Amber on Gpus. 2. Explicit Solvent Particle Mesh Ewald Journal of Chemical Theory and Computation 2013, 9, 3878-3888. (38) D.A. Case, J. T. B., R.M. Betz, D.S. Cerutti, T.E. Cheatham, III, T.A. Darden, R.E. Duke, T.J. Giese,; H. Gohlke, A. W. G., N. Homeyer, S. Izadi, P. Janowski, J. Kaus, A. Kovalenko, T.S. Lee, S. LeGrand,; P. Li, T. L., R. Luo, B. Madej, K.M. Merz, G. Monard, P. Needham, H. Nguyen, H.T. Nguyen, I.; Omelyan, A. O., D.R. Roe, A. Roitberg, R. SalomonFerrer, C.L. Simmerling, W. Smith, J. Swails,; R.C. Walker, J. W., R.M. Wolf, X. Wu, D.M. York and P.A. Kollman; University of California: San Francisco, CA, 2015. (39) Grossfield, A.; 2.0.9 ed.; University of Rochester Rochester, NY, 2013. (40) Senn, H. M.; Thiel, W. Qm/Mm Methods for Biomolecular Systems Angewandte Chemie International Edition 2009, 48, 1198-1229. (41) Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods V: Modification of Nddo Approximations and Application to 70 Elements Journal of Molecular Modeling 2007, 13, 1173-1213. (42) Walker, R. C.; Crowley, M. F.; Case, D. A. The Implementation of a Fast and Accurate Qm/Mm Potential Method in Amber Journal of Computational Chemistry 2008, 29, 1019-1031. (43) Plotnikov, N. V. Computing the Free Energy Barriers for Less by Sampling with a Coarse Reference Potential While Retaining Accuracy of the Target Fine Model J. Chem. Theory Comput. 2014, 10, 2987-3001. (44) Parkash, V.; Leppänen, V.-M.; Virtanen, H.; Jurvansuu, J. M.; Bespalov, M. M.; Sidorova, Y. A.; Runeberg-Roos, P.; Saarma, M.; Goldman, A. The Structure of the Glial Cell Line-Derived Neurotrophic Factor-Coreceptor Complex: Insights into Ret Signaling and Heparin Binding Journal of Biological Chemistry 2008, 283, 35164-35172. (45) Eigenbrot, C.; Gerber, N. X-Ray Structure of Glial Cell-Derived Neurotrophic Factor at 1.9 a Resolution and Implications for Receptor Binding Nat Struct Mol Biol 1997, 4, 435-438. Plotnikov et al. Page 38 ACS Paragon Plus Environment

Page 39 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(46) Cho, H. S.; Mason, K.; Ramyar, K. X.; Stanley, A. M.; Gabelli, S. B.; Denney, D. W., Jr.; Leahy, D. J. Structure of the Extracellular Region of Her2 Alone and in Complex with the Herceptin Fab Nature 2003, 421, 756-760. (47) Yan, B.; Steen, S.; Hambly, D.; Valliere-Douglass, J.; Bos, T. V.; Smallwood, S.; Yates, Z.; Arroll, T.; Han, Y.; Gadgil, H.et al. Succinimide Formation at Asn 55 in the Complementarity Determining Region of a Recombinant Monoclonal Antibody Igg1 Heavy Chain Journal of Pharmaceutical Sciences 2009, 98, 3509-3521. (48) Lauer, T. M.; Wood, G. P. F.; Farkas, D.; Sathish, H. A.; Samra, H. S.; Trout, B. L. Molecular Investigation of the Mechanism of Non-Enzymatic Hydrolysis of Proteins and the Predictive Algorithm for Susceptibility Biochemistry 2016, 55, 3315-3328. (49) Eyring, H. Steric Hindrance and Collision Diameters1 Journal of the American Chemical Society 1932, 54, 3191-3203. (50) Dunbrack, R. L.; Karplus, M. Backbone-Dependent Rotamer Library for Proteins Application to Side-Chain Prediction Journal of Molecular Biology 1993, 230, 543-574. (51) Leach, A. R. Molecular Modelling: Principles and Applications; Prentice Hall, 2001. (52) Luehr, N.; Ufimtsev, I. S.; Martínez, T. J. Dynamic Precision for Electron Repulsion Integral Evaluation on Graphical Processing Units (Gpus) Journal of Chemical Theory and Computation 2011, 7, 949-954. (53) Almagro, J. C.; Beavers, M. P.; Hernandez-Guzman, F.; Maier, J.; Shaulsky, J.; Butenhof, K.; Labute, P.; Thorsteinson, N.; Kelly, K.; Teplyakov, A.et al. Antibody Modeling Assessment Proteins: Structure, Function, and Bioinformatics 2011, 79, 3050-3066.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC image


Page 40 of 40

Quantifying the Risks of Asparagine Deamidation and Aspartate

Recommend Documents