==== Front Sci Rep Sci Rep Scientific Reports 2045-2322 Nature Publishing Group UK London 78136 10.1038/s41598-020-78136-1 Article Local environment effects on charged mutations for developing aggregation-resistant monoclonal antibodies Lee Jihyeon Chong Song-Ho Ham Sihyun sihyun@sookmyung.ac.kr grid.412670.60000 0001 0729 3748Department of Chemistry, The Research Institute of Natural Sciences, Sookmyung Women’s University, Cheongpa-ro 47-gil 100, Yongsan-Ku, Seoul, 04310 Korea 3 12 2020 3 12 2020 2020 10 211919 5 2020 11 11 2020 © The Author(s) 2020Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.Protein aggregation is a major concern in biotherapeutic applications of monoclonal antibodies. Introducing charged mutations is among the promising strategies to improve aggregation resistance. However, the impact of such mutations on solubilizing activity depends largely on the inserting location, whose mechanism is still not well understood. Here, we address this issue from a solvation viewpoint, and this is done by analyzing how the change in solvation free energy upon charged mutation is composed of individual contributions from constituent residues. To this end, we perform molecular dynamics simulations for a number of antibody mutants and carry out the residue-wise decomposition of the solvation free energy. We find that, in addition to the previously identified “global” principle emphasizing the key role played by the protein total net charge, a local net charge within ∼15 Å from the mutation site exerts significant effects. For example, when the net charge of an antibody is positive, the global principle states that introducing a positively charged mutation will lead to more favorable solvation. Our finding further adds that an even more optimal mutation can be done at the site around which more positively charged residues and fewer negatively charged residues are present. Such a “local” design principle accounts for the location dependence of charged mutations, and will be useful in producing aggregation-resistant antibodies. Subject terms Computational biophysicsBiophysical chemistryhttp://dx.doi.org/10.13039/501100003725National Research Foundation of KoreaNRF-2018R1D1A1B07045862NRF-2017R1A2B3010053Chong Song-Ho Ham Sihyun http://dx.doi.org/10.13039/501100014364Samsung Science and Technology FoundationSSTF-BA1401-52Ham Sihyun issue-copyright-statement© The Author(s) 2020 ==== Body Introduction Monoclonal antibodies are a growing class of biological therapeutics targeting a broad range of human disorders such as cancers and autoimmune diseases1–3. However, the development and application of antibodies as biomedicines have been severely limited due to their high propensities to aggregate under concentrated conditions4,5. Aggregation propensity is even more salient for fragments—variable heavy (VH) and light (VL) domains responsible for binding to antigens or fusions thereof—which are more preferable forms of the antibody-based regents6,7. (In the following, both the full-length antibody and the fragment shall be collectively referred to as antibody for brevity.) Several approaches have therefore been proposed to design aggregation-resistant antibodies8,9. Mutating surface residues to charged amino acids is among the promising strategies that enhances solubility and, hence, increases resistance to aggregation10–14. In this regard, the antibody’s net charge was recognized to be a key determinant of optimal charged mutations15,16: inserting positively charged residues was found to be more effective when the net charge is positive, and negatively charged residues take over the role if the net charge is negative. Since such mutations increment the magnitude of net charge, this was explained by the more strengthened electrostatic repulsion between antibody monomers, which in turn disfavors aggregation. On the other hand, the critical importance of the inserting position of charged mutations on the solubilizing activity was also recognized15–17, whose mechanism is still not well understood. This indicates the need to better understand the nature of charged mutations. The impact of charged mutations on solubilizing activity can also be rationalized from a solvation viewpoint. Indeed, the water-induced attraction is often invoked to account for biomolecular self-assembly18,19. Such an interaction can be described by the solvation free energy characterizing the affinity for the solvent water20. While it is common to apply the concept of hydrophobicity to individual amino acids, one can argue the “overall protein hydrophobicity” in terms of the solvation free energy defined for a whole protein. In fact, we have recently demonstrated that the overall protein hydrophobicity is the key factor controlling protein aggregation propensity21. The same “global” design principle mentioned above—introducing positively (negatively) charged mutations enhances the aggregation resistance when the net charge is positive (negative)—has also been derived from the observation that the solvent water strikingly discriminates positive and negative charges depending on the protein net charge. This solvation perspective was recently adopted by Schäfer and coworkers22, who have successfully designed antibody mutants of improved solubility guided by the solvation thermodynamics analysis. On the other hand, it was also found that the variation in solvation free energy upon charged mutation markedly depends on the insertion location, which is outside the scope of the global principle. An additional guiding principle that takes into account the spatial neighbors of the mutation site is therefore necessary. In this paper, we explore such a “local” design principle for charged mutations from detailed solvation thermodynamics analyses. We investigate the same antibody systems studied in Ref.22 since improved solution-state properties of those systems through the optimization of the solvation free energy have been experimentally confirmed. We also analyze the single-domain antibody for which the solubilizing activity of charged mutations was experimentally measured16. We conduct molecular dynamics simulations and carry out solvation free energy calculations for these systems. We then perform the residue-wise decomposition analysis of the solvation free energy23–25. This allows us to gain detailed knowledge on how the solvation free energy change is composed of individual contributions from constituent residues. We particularly focus on the contributions from the spatial neighbors around the mutation site, whose understanding will provide us with new insights on the inserting-location dependence of charged mutations. Thereby, we would like to derive a local design principle that will be useful in producing aggregation-resistant antibodies. Methods System description We investigated the Fv fragments consisting of VH and VL domains of two monoclonal antibodies (mAbs), to be referred to as mAb1 and mAb2 (Fig. 1a,b). The former (mAb1) is the antibody that neutralizes human immunodeficiency virus-type 1 (HIV-1)26, and the latter (mAb2) is a human anti-DNA autoantibody27. The starting structure of mAb1 was taken from the X-ray study (PDB entry 3RU8)28. For mAb2, we built a homology model since no experimental structure is available. This was done with the antibody modeler module in the Molecular Operating Environment (MOE) software29 using the structure of PDB entry 1DFB30 as a template. We analyzed 10 mutants of mAb1 and 8 mutants of mAb2 introduced in Ref.22 (see Tables 1 and 2). The mutation sites were chosen from the solvent-exposed residues by excluding conserved residues and those within or near the complementarity-determining regions (CDRs; see Fig. 1a,b).Figure 1 (a, b) Structures of the Fv fragments comprising variable heavy (VH) and light (VL) domains of (a) the wild-type mAb1 (PDB entry 3RU8) and (b) the wild-type mAb2 (modeled based on PDB entry 1DFB). The mutation sites studied in the present work are indicated by stick representations and green labels. Heavy and light complementarity-determining regions (CDRs) are colored blue and red, respectively. (c) Structure of the single-domain (VH) antibody Aβ18-27 dAb (modeled based on PDB entry 3B9V). In the “wild type” dAb, three Ala residues (stick representations colored green) are substituted to the mutation sites (enclosed by a green dashed oval) located at the N-terminus of the third CDR (CDR3). PyMOL version 1.8.2 (https://pymol.org) was used to generate protein figures. Table 1 Solvation free energy changes upon mutating mAb1. ΔQa ΔGsolv (kcal/mol)b ΔGsolv from Ref.22 Heavy-chain mutations  E10G +1 -128.2 ± 0.9 -135.6 ± 20.5  D73N +1 -85.1 ± 16.7 -87.4 ± 19.3  A76K +1 -121.8 ± 14.2 -116.1 ± 21.1  E10G/D73N/A76K +3 -343.8 ± 8.4 -396.9 ± 20.1 Light-chain mutations D60S +1 -73.9 ± 15.3 -102.7 ± 21.5 E80Q +1 -110.4 ± 4.6 -106.9 ± 22.5 D60S/E80Q +2 -160.9 ± 9.8 -223.7 ± 21.6 Heavy- and light-chain mutations E10G/A76K/E80Q +3 -268.0 ± 13.3 -285.7 ± 21.5 D73N/A76K/E80Q +3 -282.4 ± 9.8 -292.3 ± 22.0 E10G/D73N/A76K/D60S/E80Q +5 -503.1 ± 3.5 -484.0 ± 21.1 aIncrement in charge upon mutation(s); bΔGsolv=Gsolv(mutant)-Gsolv(wild type). Table 2 Solvation free energy changes upon mutating mAb2. ΔQa ΔGsolv [kcal/mol]b ΔGsolv from Ref.22 Heavy-chain mutations Q13K +1 +12.0 ± 13.3 -1.5 ± 13.5 D73N +1 -37.0 ± 7.0 -82.2 ± 13.1 Q115K +1 -43.6 ± 12.0 -40.2 ± 13.1 Q13K/D73N/Q115K +3 -218.1 ± 11.4 -216.9 ± 14.1 Light-chain mutations D32Y +1 +60.2 ± 19.5 +56.3 ± 12.5 S77R +1 +22.8 ± 11.9 -31.6 ± 13.8 D32Y/S77R +2 -12.2 ± 7.1 +10.1 ± 13.3 Heavy- and light-chain mutations Q13K/D73N/Q115K/D32Y/S77R +5 -316.1 ± 20.7 -205.6 ± 13.7 aIncrement in charge upon mutation(s); bΔGsolv=Gsolv(mutant)-Gsolv(wild type). Table 3 Solvation free energy changes upon mutating Aβ18-27 dAb. ΔQa ΔGsolv [kcal/mol]b Aggregationc RRR-Aβ18-27 +3 -8.5 ± 4.5 Present DDD-Aβ18-27 -3 -507.7 ± 4.0 Eliminated aIncrement in charge upon mutation(s); bΔGsolv=Gsolv(mutant)-Gsolv(wild type); cPresence/absence of the aggregation as reported in Ref.16. We also studied the single-domain (VH) antibody (dAb), to be referred to as Aβ18-27 dAb, which incorporates hydrophobic residues 18-VFFAEDVGSN-27 taken from the Alzheimer’s amyloid-β (Aβ) peptide within the third CDR (CDR3; see Fig. 1c). Aβ18-27 dAb binds to Aβ oligomers and fibrils with submicromolar affinity31, but it is also prone to self-aggregate within days at 25 ∘C16. We analyzed 2 mutants of Aβ18-27 dAb in which three charged residues are inserted at the N-terminus of CDR3, whose solubilizing activity was measured experimentally (Table 3)16. The starting structure of Aβ18-27 dAb was also modeled with the MOE software29 using the structure of PDB entry 3B9V32 as a template. Molecular dynamics simulations We used the pmemd.cuda module in AMBER18 package33 to perform explicit-water molecular dynamics simulations for all the systems. The ff14SB force field34 and TIP3P model35 were employed for proteins and water, respectively. The simulations were done under neutral pH, where Glu and Asp carry a negative charge and Lys and Arg a positive charge. The resulting net charges of mAb1, mAb2 and Aβ18-27 dAb are + 6, + 3 and 0, respectively. The whole charge of each simulation system was neutralized by counter Cl- ions, and additional Na+ and Cl- were included to achieve a 150 mM ionic concentration. We applied the particle mesh Ewald method36 to handle long-range Coulomb interactions, and short-range interactions were treated by a 10 Å cutoff. Berendsen’s thermostat and barostat37 were used for constant temperature and pressure (300 K and 1 bar). Two independent 100 ns production simulations were carried out for each system. Solvation free energy analysis Based on the general expression for the solvation free energy Gsolv known as the Kirkwood charing formula, we have derived the following exact atomic decomposition of Gsolv23: 1 Gsolv=∑αGsolv,α with 2 Gsolv,α=4π∑γ∫01dλ∫r2dr∂uαγ(r;λ)∂λργgαγ(r;λ) Here α and γ refer to the solute (protein) and solvent sites, respectively; λ is the coupling parameter that introduces the solute-solvent interaction uαγ(r) such that uαγ(r;λ=0)=0 and uαγ(r;λ=1)=uαγ(r); ργ is the average number density of site γ; and gαγ(r;λ) is the solute-solvent radial distribution function corresponding to uαγ(r;λ). (In the original work23, λ consists of two parameters, λ=(λ1,λ2), to separately control the short-range and Coulomb interactions, but such a complexity is suppressed here for simplicity.) By an appropriate grouping of contributions from constituent atoms (Gsolv,α), one obtains a residue-wise decomposition of Gsolv. In the present work, the protein-solvent distribution function was computed from the three-dimensional reference interaction site model (3D-RISM) theory38,39 (see the Supplementary Information for details). Results Structural and solvation free energy changes upon mutations in mAb1 and mAb2 We performed explicit-water molecular dynamics simulations for the wild-type and 10 mutants of mAb1 (Table 1) and for the wile-type and 8 mutants of mAb2 (Table 2). Protein structures were stable during the simulation time (100 ns) in all the systems: the Cα RMSD (root-mean-square deviation) from the respective initial structure stayed within ∼3 Å. We then computed the solvation free energy Gsolv for each system using the simulated structures sampled with a 1 ns interval. In the following, we will focus on the change in solvation free energy upon mutation defined by ΔGsolv=Gsolv(mutant)-Gsolv(wildtype). The average and standard error of ΔGsolv were estimated based on the two independent production runs. The results for mAb1 and mAb2 are presented in Tables 1 and 2, along with the comparison with the previous work. (Results of individual trajectories are reported in Supplementary Table S1 and S2.) Overall, our numerical results for ΔGsolv are very similar to those reported previously22, though some differences are discernible. In particular, ΔGsolv values for the Q13K, S77R and D32Y/S77R mutants of mAb2 show the opposite signs compared to the previous work (Table 2). To examine such numerical differences and how they affect the main points (in particular, a local design principle) of the present work, we have carried out two additional independent 100 ns simulations for these systems, and the results are summarized in Supplementary Tables S3 and S4. As demonstrated there, the results from the additional simulations are in better agreement with the previous work, especially concerning the signs of ΔGsolv values. This indicates that, with two independent 100 ns simulations, there still remains non-negligible uncertainty in numerical results for ΔGsolv. We will come back to this issue later in Discussion section, where we argue that our main points are nevertheless not significantly altered. Single-point mutations in mAb1 We start from analyzing the single-point mutants of mAb1. Individual residue contributions to ΔGsolv upon the E10G (heavy chain) mutation are shown in Fig. 2a; those upon the A76K (heavy chain) mutation in Fig. 2b; and the results for the D73N (heavy chain), D60S (light chain) and E80G (light chain) mutations are displayed in Supplementary Fig. S1. The single-point mutations considered here are either the mutation of a negatively charged residue to a neutral one or the mutation of a neutral residue to a positively charged one, which increments the system charge by +1. (We recall here that the total charge of the wild-type mAb1 is + 6.) Hence, according to the aforementioned global principle21, positively charged residues get more favorably solvated (resulting in more negative ΔGsolv), whereas negatively charged residues become less favorably solvated (more positive ΔGsolv). This explains why the contributions to ΔGsolv from the positively charged residues (colored blue) tend to exhibit negative changes and those from the negatively charged residues (colored red) show the opposite trend.Figure 2 (a, b) Residue-wise decomposition of ΔGsolv for (a) E10G and (b) A76K mutants of mAb1 (blue, red and black colors refer to positively charged, negatively charged and neutral residues, respectively). The mutation site is indicated by the green arrow. dm (in Å) denotes the distance to the mutation site as illustrated in the inset (green stick representation, mutation site; cyan and orange representations, positively and negatively charged residues, respectively). Residues involved in the formation/breaking of salt-bridges upon mutation are represented by SB +/- (green if the mutation site is involved, and black otherwise). PyMOL version 1.8.2 (https://pymol.org) was used to generate protein figures. On the other hand, the presence of charged residues that exhibit more pronounced variations is discernible (indicated by the dashed arrows or SB +/- in Fig. 2 and Supplementary Fig. S1). We first recognize that such large variations tend to occur in the same chain into which the mutation is introduced: for example, in the heavy-chain E10G mutation (Fig. 2a), most of the large changes show up in the heavy chain, but not in the light chain. This indicates that the distance to the mutation site matters. We therefore added the distance information (dm in Å; the minimum side-chain heavy atom distance to the mutation site) in Fig. 2 and Supplementary Fig. S1. In fact, we find that the residues exhibiting the large solvation free energy charges are located within ∼15 Å from the mutation site. (In this regard, we notice that those distances, and the insets in Fig. 2, were obtained using the wile-type mAb1 structure. This is because we are interested in recipes for designing aggregation-resistant mutants that can be utilized solely based on the knowledge of the wild-type protein.) Second, we observe a sizable difference in the solvation free energy change at the mutation site (indicated by the green arrow) between the E10G (Fig. 2a) and A76K (Fig. 2b) mutants. This implies that the contents of the mutation (such as a deletion of negatively charged residue versus an addition of positively charged residue) influence ΔGsolv value at the mutation site. Finally, it is seen that the residues involved in the formation/breaking of salt-bridges (indicated by SB +/-, respectively) upon mutation are significantly affected. However, except for those directly associated with the mutation site (SB +/- colored green), it is difficult to predict in advance the formation/breaking of salt-bridges far from the mutation site (SB +/- colored black). Therefore, those salt-bridges will not be taken into consideration in our exploration of design principles. The key to understanding the pronounced variations in ΔGsolv near and at the mutation site lies in the solvation structure21. Since the net charge of mAb1 is positive (+ 6), the protein-water electrostatic interaction induces such a long-distance orientational distribution of water where the dipole moment is directed outward from the protein. Those water molecules in turn produce an electrostatic potential that energetically favors positive charges and disfavors negative ones on the protein surface. When a charged mutation of +1 is introduced, such an electrostatic potential strengthens around the mutation site. This explains the large variations in ΔGsolv, including their signs, of the residues located within ∼15 Å from the mutation site. When this mutation is a deletion of a negatively charged residue, there is no charged residue that is affected by the strengthened electrostatic potential at the mutation site. On the other hand, if the mutation is an addition of a positively charged residue, the mutation site fully appreciates the strengthened electrostatic potential. This accounts for the large difference in ΔGsolv values at the mutation site between the E10G and A76K mutants. The results for the other mutants shown in Supplementary Fig. S1 can be rationalized in a similar manner. Single-point mutations in mAb2 We next investigate the single-point mutants of mAb2. In contrast to large negative ΔGsolv values (-73.9 to -128.2 kcal/mol; see Table 1) exhibited by the single-point mutants of mAb1, corresponding ΔGsolv values of mAb2 are much less negative (the most negative value is -43.6 kcal/mol) and even take positive values (Table 2). The less negative ΔGsolv for mAb2 can be understood from the smaller net charge of the wild-type mAb2 (+3) than that of mAb1 (+6). Indeed, according to the simple continuum Born model40, the solvation free energy change caused by increasing a charge from Q to Q+ΔQ (ΔQ>0) is proportional to -(2QΔQ+ΔQ2), and it is negative when Q>0 and it becomes more negative for larger Q. Thus, the positive solvation free energy change upon increasing the net charge observed for the mutants of mAb2 is counter-intuitive. To understand such a nontrivial solvation behavior of mAb2, we analyze in detail the Q13K (heavy chain) and Q115K (heavy chain) mutants. These are the same type of mutation (glutamine to lysine), but exhibit strikingly contrasting ΔGsolv values (+12.0 and -43.6 kcal/mol, respectively) depending on the mutation location, and, hence, constitute good examples for elucidating the position dependence of charged mutations. Residue-decomposed ΔGsolv values for these mutants are shown in Fig. 3a,b. Main features of Fig. 3a,b can be rationalized in the same way as we explained above for mAb1: since mAb2’s net charge is positive, positive (negative) residues tend to exhibit negative (positive) solvation free energy changes; more pronounced effects show up in those residues located within ∼15 Å from the mutation site; and the ΔGsolv value at the mutation site is large negative since a new positively charged residue is introduced there. However, an inspection of Fig. 3a for Q13K and Fig. 3b for Q115K reveals a difference in the distribution of charged residues around (<∼15 Å) the respective mutation site: the Q13K mutation site is surrounded by two negatively charged residues, whereas the Q115K mutation site by two negatively charged residues and three positively charged residues. Those nearby negatively charged residues provide significant positive contributions to ΔGsolv, and this accounts for the positive ΔGsolv of the Q13K mutant. On the other hand, in the case of Q115K, such large positive contributions to ΔGsolv arising from the nearby negatively charged residues are sufficiently compensated by large negative contributions from the nearby positively charged residues. Thus, a local net charge within ∼15 Å from the mutation site exerts significant impacts on the ΔGsolv value. This also provides an explanation of the location dependence of charged mutation from the solvation perspective.Figure 3 (a, b) Residue-wise decomposition of ΔGsolv for (a) Q13K and (b) Q115K mutants of mAb2 (blue, red and black colors refer to positively charged, negatively charged and neutral residues, respectively). The mutation site is indicated by the green arrow. dm (in Å) denotes the distance to the mutation site as illustrated in the inset (green stick representation, mutation site; cyan and orange representations, positively and negatively charged residues, respectively). Residues involved in the formation/breaking of salt-bridges upon mutation are represented by SB +/- (green if the mutation site is involved, and black otherwise). PyMOL version 1.8.2 (https://pymol.org) was used to generate protein figures. Individual residue contributions to ΔGsolv of the other single-point mutants of mAb2—D73N (heavy chain; ΔGsolv=-37.0 kcal/mol), D32Y (light chain; ΔGsolv=+60.2 kcal/mol), and S77R (light chain; ΔGsolv=+22.8 kcal/mol) shown in Supplementary Fig. S2—can be rationalized in a similar manner. ΔGsolv for the D73N mutant is largely negative since the mutation site is surrounded by three positively charged residues. ΔGsolv for the D32Y mutant is positive mainly because a bulky side chain is introduced at the solvent exposed residue, which produces a large positive ΔGsolv at the mutation site. This suggests that such a mutation that introduces a solvent-exposed bulky side chain should be avoided. ΔGsolv for the S77R mutant is positive since the mutation site is surrounded by two positively charged residues and three negatively charged residues. Multipoint mutations in mAb1 and mAb2 Finally, we deal with the multipoint mutations. Because of the larger (≥+2) increment of the net charge in these mutations, the protein-water electrostatic effects on ΔGsolv become stronger. Therefore, positive (negative) residues exhibit more pronounced negative (positive) solvation free energy changes than those observed in the single-point mutations. This is seen, e.g., in Fig. 4a showing individual residue contributions to ΔGsolv of the E10G/D73N/A76K (heavy chain) mutant of mAb1. We also notice from this figure that the solvation free energy changes induced by the three mutations can roughly be described as a superposition of the effects caused by individual mutations: the five residues strongly affected by the single E10G mutation (indicated by the dashed arrows in Fig. 2a) are shown here with the cyan arrows; the three residues associated with the A76K mutation (taken from Fig. 2b) are represented by the magenta arrows; and the three residues whose variations are large in the D73N mutant (Supplementary Fig. S1) are indicated by the orange arrows. (A certain complexity arises here since a residue that exhibits a large solvation free energy change upon one mutation may be deleted in another mutation. Such a case is exemplified by the dashed magenta arrow in Fig. 4a.) The corresponding result for the D60S/E80Q mutations in the light chain of mAb1 (Fig. 4b) can be rationalized likewise in terms of the individual D60S and E80Q mutations shown in Supplementary Fig. S1. The result for the E10G/D73N/A76K/D60S/E80Q mutations (Fig. 4c) can also be roughly described as a superposition of those shown in Fig. 4a,b, and this applies to the other (E10G/A76K/E80Q and D73N/A76K/E80Q) multipoint mutants of mAb1 (Supplementary Fig. S3).Figure 4 (a–c) Residue-wise decomposition of ΔGsolv for (a) E10G/D73N/A76K, (b) D60S/E80Q, and (c) E10G/D73N/A76K/D60S/E80Q mutants of mAb1 (blue, red and black colors refer to positively charged, negatively charged and neutral residues, respectively). The mutation sites are indicated by the green arrows. Charged residues that exhibit pronounced variations in individual single-point mutations (taken from Fig. 2 and Supplementary Fig. S1) are represented by the solid arrows labeled with respective mutations. When those charged residues are deleted upon multipoint mutations, the solid arrows are replaced by the dashed ones. The residue-wise decomposition of ΔGsolv for the multipoint mutants of mAb2, shown in Fig. 5, can be understood similarly. The solvation free energy change for the D32Y/S77R mutant (ΔGsolv=-12.2 kcal/mol) is still not large negative in spite of the double-point charged mutation. This is not only because the bulky side chain (from D32Y) is involved, but also because five positively charged residues and six negatively charged residues are present close to those mutation sites (Fig. 5b). Upon introducing additional charged mutations, more positively charged residues from the Q13K/D73N/Q115K mutations (Fig. 5a) come into play, and the solvation free energy change for the Q13K/D73N/Q115K/D32Y/S77R mutant becomes a large negative value (Fig. 5c).Figure 5 (a–c) Residue-wise decomposition of ΔGsolv for (a) Q13K/D73N/Q115K, (b) D32Y/S77R, and (c) Q13K/D73N/Q115K/D32Y/S77R mutants of mAb2 (blue, red and black colors refer to positively charged, negatively charged and neutral residues, respectively). The mutation sites are indicated by the green arrows. Charged residues that exhibit pronounced variations in individual single-point mutations (taken from Fig. 3 and Supplementary Fig. S2) are represented by the solid arrows labeled with respective mutations. Multipoint mutations in Aβ18-27 dAb Now let us turn our attention to Aβ18-27 dAb. We studied the RRR-Aβ18-27 and DDD-Aβ18-27 mutants in which three charged residues (RRR or DDD) are inserted at the N-terminus of CDR3 (Fig. 1c). For the computation of ΔGsolv and its residue-wise decomposition, we also analyzed AAA-Aβ18-27 (i.e., three alanine residues are substituted to the mutation sites) which was considered as the wild-type. We carried out two independent 100 ns molecular dynamics simulations for these systems, and the simulation structures were stable for all the systems (Cα-RMSD stayed within ∼3 Å). We then computed the solvation free energy change ΔGsolv=Gsolv(mutant)-Gsolv(wildtype), and the results are summarized in Table 3 (see also Supplementary Table S5). The net charge of the wild-type Aβ18-27 dAb is 0, and its change is symmetrical for RRR-Aβ18-27 (+3) and DDD-Aβ18-27 (-3). Therefore, the global principle21 cannot tell us which of the charged mutations leads to more favorable solvation. (This is the primary reason why this particular dAb was chosen among several related dAbs studied in Refs.15,16.) Interestingly, the solubilizing activity of RRR and DDD mutations was found to be contrastive; experimental measurements indicate that aggregation is still observable for RRR-Aβ18-27, but it is eliminated for DDD-Aβ18-2716. Correspondingly, we find that ΔGsolv for DDD-Aβ18-27 is strikingly more negative than the one for RRR-Aβ18-27 (Table 3). Thus, there must be some local environment effects that distinguish positively and negatively charged residues at the N-terminus of CDR3. The residue-wise decomposition analysis of ΔGsolv is suited to address such local effects, and the results for RRR-Aβ18-27 and DDD-Aβ18-27 are shown in Fig. 6. We first note that, since the total charge of RRR-Aβ18-27 is positive, positive and negative residues respectively tend to exhibit negative and positive ΔGsolv values (Fig. 6a), whereas this trend is inverted for DDD-Aβ18-27 (Fig. 6b). We observe from Fig. 6a,b and the insets that more negatively charged residues are present than positively charged ones within ∼15 Å from the mutation sites. This explains why the ΔGsolv values at the mutation sites in RRR-Aβ18-27 is much less negative than those in DDD-Aβ18-27. Thus, it is the local principle, derived above from the analyses of mAb1 and mAb2, that accounts for the distinct solubilizing activity of RRR and DDD mutations.Figure 6 (a, b) Residue-wise decomposition of ΔGsolv for (a) RRR-Aβ18-27 and (b) DDD-Aβ18-27 mutants of Aβ18-27 dAb (blue, red and black colors refer to positively charged, negatively charged and neutral residues, respectively). The mutation sites are indicated by the green arrow. dm (in Å) denotes the distance to the mutation sites as illustrated in the inset (green dashed oval, mutation sites; cyan and orange stick representations, positively and negatively charged residues, respectively). Residues involved in the formation/breaking of salt-bridges upon mutation are represented by SB +/- (green if the mutation sites are involved, and black otherwise). PyMOL version 1.8.2 (https://pymol.org) was used to generate protein figures. Discussion Protein aggregation has been an area of major focus in its own right because of its relation to human diseases41. Factors promoting aggregation-prone nature have therefore been intensively investigated, and several algorithms have come out to rationalize and predict protein aggregation propensity42–45. These algorithms mainly focus on sequence characteristics such as the β-sheet forming propensity and the amino acid hydrophobicity. They have also been applied to identify aggregation-prone regions in antibodies, and it was found that such regions are typically located within CDRs responsible for antigen binding46–48. On the other hand, we have recently proposed a different type of strategy from a solvation viewpoint21. This is natural since whether a protein becomes aggregation-prone or remains soluble in aqueous environments should critically depend on the affinity for water. In fact, it was demonstrated that the solvation free energy is the key factor controlling the protein aggregation propensity21. Since the solvation free energy is largely affected by charged residues, we could thereby derive a design principle for increasing the resistance to aggregation through charged mutations. However, the previously derived principle is a global one in the sense that it takes into consideration only the global net charge of a protein and the sign (positive or negative) of charged mutations. It therefore cannot discriminate, e.g., the Q13K and Q115K mutants of mAb2 studied in the present work, which however exhibit strikingly contrasting solvation free energy changes (+12.0 and -43.6 kcal/mol, respectively; see Table 2). Through the residue-wise decomposition analysis of the solvation free energy changes, it is elucidated that the contrasting behavior of the Q13K and Q115K mutants is caused by the difference in local charge surrounding the mutation site. Thus, the spatial distribution of nearby (<∼ 15 Å) charged residues is found to significantly impact the solvation free energy change. We also observe that the effects of multipoint mutations can be roughly described as a superposition of those from single-point mutations. Furthermore, our analysis for mAb2 reveals why multiple mutations (more than two) are sometimes necessary to confer sufficiently negative solvation free energy. These results are in accord with the experimental observation that multiple mutations are necessary to considerably improve aggregation resistance and their effects are largely additive17. Antibodies of neutral net charge are also the systems for which the previously derived global principle alone is insufficient: it cannot propose which type (positive or negative) of charged mutations should be introduced to enhance the solubility. In this connection, it was observed experimentally that inserting negatively charged residues (DDD) is more beneficial than inserting positively charged residues (RRR) in preventing the aggregation of Aβ18-27 dAb of neutral net charge16. This observation cannot be rationalized by the electrostatic repulsion since its magnitude is the same between positive charges and between negative charges. On the other hand, we demonstrate that the solvation free energy clearly discriminates DDD-Aβ18-27 and RRR-Aβ18-27 and that it can be accounted for by the local environment surrounding the mutation sites. This indicates a potential advantage of adopting the solvation viewpoint in arguing protein aggregation. Finally, let us make comments related to the convergence of solvation free energy calculations. As stated at the beginning of Results section, additional two independent 100 ns simulations were carried out for the Q13K, S77R and D32Y/S77R mutants of mAb2, and non-negligible differences were found between the average ΔGsolv values computed from the original and additional simulations (Supplementary Table S3). To examine the discrepancies, we plotted ΔGsolv as a function of simulation time for individual trajectories (Supplementary Figs. S4, S5 and S6). We observe that ΔGsolv values fluctuate significantly (the standard deviation in individual trajectories is ∼70 to 90 kcal/mol). Such large fluctuations in individual trajectories reflect the fact that the solvation free energy depends on fine details of the underlying protein structure: for example, the formation of just a single salt-bridge in a protein can lead to ∼50 kcal/mol variation23. Therefore, even subtle differences in the population of hydrogen bonds and salt-bridges in different trajectories can significantly affect the average ΔGsolv value, and it is extremely difficult to obtain its converged result (e.g., within an error of ∼1 kcal/mol). The average ΔGsolv values reported in the present work were computed from relatively short (100 ns) simulations, with which typical errors are of the order of ∼10 kcal/mol (Tables 1, 2 and 3). Hence, the average ΔGsolv values for the Q13K, S77R and D32Y/S77R mutants of mAb2 are not decisive from our calculations: ΔGsolv∼±10 kcal/mol for these mutants with errors of the order of ∼10 kcal/mol. Nevertheless, this observation does not invalidate the present work. First, the global and local effects on ΔGsolv upon charged mutation—the main points of our work—are not affected. This is also demonstrated in Supplementary Figs. S4, S5 and S6, where the residue-wise decompositions of ΔGsolv computed from the original and additional simulations are compared. We observe that the global effect is conserved: positive (negative) residues tend to exhibit negative (positive) solvation free energy changes in both the the original and additional simulations, and this can be explained by the fact that the net charge of the mutants under consideration is positive. In addition, the same local effects are observed in both the the original and additional simulations: the charged residues within ∼15 Å from the mutation site(s) are affected in the same manner. Second, our goal is to provide design principles of aggregation-resistant antibodies upon charged mutations. As stated above, it is already recognized that multiple charged mutations are necessary to improve the aggregation resistance17. For the multiple charged mutations, the typical magnitude of ΔGsolv significantly exceeds 100 kcal/mol (see Tables 1, 2 and 3), which is reliable even with errors of ∼10 kcal/mol. Thus, although there is a problem if our work is concerned with the precise prediction of ΔGsolv, the convergence of the global and local effects is attained already with the simulations reported here and our computational methods are reliably applicable in designing aggregation-resistant antibodies. (Mutants for which ΔGsolv∼±10 kcal/mol, such as Q13K, S77R and D32Y/S77R mentioned above, can simply be excluded from candidates of aggregation-resistant antibodies). Conclusions In summary, we explore here strategies for designing aggregation-resistant antibodies by charged mutations. This is done on the basis of a solvation thermodynamics perspective rather than focusing on the sequence characteristics. By combining our previous work and the analyses presented here, we identify three critical factors that significantly impact charged mutations: (1) the total net charge of an antibody; (2) the contents of mutation (the choice of the positively or negatively charged residues and avoiding the use of residues of bulky side chains); and (3) the local net charge within ∼15 Å from the mutation site. The third factor accounts for the location dependence of charged mutations, and needs to be taken into consideration in searching for an optimal insertion location. A prominent advantage of charged mutations is that they can be inserted outside the complementarity-determining regions mediating target binding, i.e., the antibody solubility can in principle be improved without reducing the binding affinity. We hope that our design principles will be useful in producing aggregation-resistant antibody therapeutics. Supplementary information Supplementary Information Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Supplementary information is available for this paper at 10.1038/s41598-020-78136-1. Acknowledgements This work was supported by the Samsung Science and Technology Foundation under Project Number SSTF-BA1401-52 and the National Research Foundation of Korea (NRF) (NRF-2017R1A2B3010053 and NRF-2018R1D1A1B07045862). Author contributions J.L., S.-H.C. and S.H. conducted the research and wrote the manuscript. Competing interests The authors declare no competing interests. ==== Refs References 1. Nelson AL Dhimolea E Reichert JM Development trends for human monoclonal antibody therapeutics Nat. Rev. Drug Discov. 2010 9 767 774 10.1038/nrd3229 20811384 2. Sliwkowski MX Mellman I Antibody therapeutics in cancer Science 2013 341 1192 1198 10.1126/science.1241145 24031011 3. Carter PJ Lazar GA Next generation antibody drugs: Pursuit of the high-hanging fruit Nat. Rev. Drug Discov. 2018 17 197 223 10.1038/nrd.2017.227 29192287 4. Lowe D Aggregation, stability, and formulation of human antibody therapeutics Adv. Protein Chem. Struct. Biol. 2011 84 41 61 10.1016/B978-0-12-386483-3.00004-5 21846562 5. Li W Antibody aggregation: Insights from sequence and structure Antibodies 2016 5 19 10.3390/antib5030019 6. Jespers L Schon O Famm K Winter G Aggregation-resistant domain antibodies selected on phage by heat denaturation Nat. Biotechnol. 2004 22 1161 1165 10.1038/nbt1000 15300256 7. Holliger P Hudson PJ Engineered antibody fragments and the rise of single domains Nat. Biotechnol. 2005 23 1126 1136 10.1038/nbt1142 16151406 8. Lee CC Perchiacca JM Tessier PM Toward aggregation-resistant antibodies by design Trends Biotechnol. 2013 31 612 620 10.1016/j.tibtech.2013.07.002 23932102 9. Roberts CJ Therapeutic protein aggregation: Mechanisms, design and control Trends Biotechnol. 2014 32 372 380 10.1016/j.tibtech.2014.05.005 24908382 10. Lawrence MS Phillips KJ Liu DR Supercharging proteins can impart unusual resilience J. Am. Chem. Soc. 2007 129 10110 10112 10.1021/ja071641y 17665911 11. Wu S-J Structure-based engineering of a monoclonal antibody for improved solubility Protein Eng. Des. Sel. 2010 23 643 651 10.1093/protein/gzq037 20543007 12. Simeonov P Berger-Hoffmann R Hoffmann R Sträter N Zuchner T Surface supercharged human enteropeptidase light chain shows improved solubility and refolding yield Protein Eng. Des. Sel. 2011 24 261 268 10.1093/protein/gzq104 21084283 13. Perchiacca JM Bhattacharya M Tessier PM Mutational analysis of domain antibodies reveals aggregation hotspots within and near the complementarity determining regions Proteins 2011 79 2637 2647 10.1002/prot.23085 21732420 14. Miklos AE Structure-based design of supercharged, highly thermoresistant antibodies Chem. Biol. 2012 19 449 455 10.1016/j.chembiol.2012.01.018 22520751 15. Perchiacca JM Ladiwala ARA Bhattacharya M Tessier PM Aggregation-resistant domain antibodies engineered with charged mutations near the edges of the complementarity-determining regions Protein Eng. Des. Sel. 2012 25 591 601 10.1093/protein/gzs042 22843678 16. Perchiacca JM Lee CC Tessier PM Optimal charged mutations in the complementarity-determining regions that prevent domain antibody aggregation are dependent on the antibody scaffold Protein Eng. Des. Sel. 2014 27 29 39 10.1093/protein/gzt058 24398633 17. Dudgeon K General strategy for the generation of human antibody variable domains with increased aggregation resistance Proc. Natl. Acad. Sci. USA 2012 109 10879 10884 10.1073/pnas.1202866109 22745168 18. Kauzmann W Some factors in the interpretation of protein denaturation Adv. Protein. Chem. 1959 14 1 64 10.1016/S0065-3233(08)60608-7 14404936 19. Tanford C The hydrophobic effect and the organization of living matter Science 1978 200 1012 1018 10.1126/science.653353 653353 20. Ben-Naim A Hydrophobic Interactions 1980 New York Plenum Press 21. Chong S-H Ham S Interaction with the surrounding water plays a key role in determining the aggregation propensity of proteins Angew. Chem. Int. Ed. 2014 53 3961 3964 10.1002/anie.201309317 22. Kuhn AB Improved solution-state properties of monoclonal antibodies by targeted mutations J. Phys. Chem. B 2017 121 10818 10827 10.1021/acs.jpcb.7b09126 29135256 23. Chong S-H Ham S Atomic decomposition of the protein solvation free energy and its application to amyloid-beta protein in water J. Chem. Phys. 2011 135 034506 10.1063/1.3610550 21787012 24. Chong S-H Ham S Site-directed analysis on protein hydrophobicity J. Comput. Chem. 2014 35 1364 1370 10.1002/jcc.23631 24817476 25. Chong S-H Ham S Distinct role of hydration water in protein misfolding and aggregation revealed by fluctuating thermodynamics analysis Acc. Chem. Res. 2015 48 956 965 10.1021/acs.accounts.5b00032 25844814 26. Burton DR Efficient neutralization of primary isolates of HIV-1 by a recombinant human monoclonal antibody Science 1994 266 1024 1027 10.1126/science.7973652 7973652 27. Rioux JD Zdárský E Newkirk MM Rauch J Anti-DNA and anti-platelet specificities of SLE-derived autoantibodies: Evidence for CDR2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{{\rm h}}$$\end{document} h mutations and CDR3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{{\rm h}}$$\end{document} h motifs Mol. Immunol. 1995 32 683 696 10.1016/0161-5890(95)00045-G 7659095 28. Azoitei ML Computation-guided backbone grafting of a discontinuous motif onto a protein scaffold Science 2011 334 373 376 10.1126/science.1209368 22021856 29. Maier JK Labute P Assessment of fully automated antibody homology modeling protocols in molecular operating environment Proteins 2014 82 1599 1610 10.1002/prot.24576 24715627 30. He XM Rüker F Casale E Carter DC Structure of a human monoclonal antibody Fab fragment against gp41 of human immunodeficiency virus type 1 Proc. Natl. Acad. Sci. USA 1992 89 7154 7158 10.1073/pnas.89.15.7154 1496010 31. Perchiacca JM Ladiwala ARA Bhattacharya M Tessier PM Structure-based design of conformation- and sequence-specific antibodies against amyloid \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} β Proc. Natl. Acad. Sci. USA 2012 109 84 89 10.1073/pnas.1111232108 22171009 32. Barthelemy PA Comprehensive analysis of the factors contributing to the stability and solubility of autonomous human VH domains J. Biol. Chem. 2008 283 3639 3654 10.1074/jbc.M708536200 18045863 33. Case DA AMBER 18 2018 San Francisco University of California 34. Maier JA ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB J. Chem. Theory Comput. 2015 11 3696 3713 10.1021/acs.jctc.5b00255 26574453 35. Jorgensen WL Chandrasekhar J Madura JD Impey RW Klein ML Comparison of simple potential functions for simulating liquid water J. Chem. Phys. 1983 79 926 935 10.1063/1.445869 36. Darden T York D Pedersen L Particle mesh Ewald: An \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n$$\end{document} n -\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log (n)$$\end{document} log ( n ) method for Ewald sums in large systems J. Chem. Phys. 1993 98 10089 10092 10.1063/1.464397 37. Berendsen HJC Postma JPM van Gunsteren WF DiNola A Haak JR Molecular dynamics with coupling to an external bath J. Chem. Phys. 1984 81 3684 3690 10.1063/1.448118 38. Kovalenko A Hirata F Three-dimensional RISM theory for molecular liquids and solid-solid interfaces Molecular Theory of Solvation 2003 Dordrecht Kluwer Academic 169 275 39. Imai T Harano Y Kinoshita M Kovalenko A Hirata F A theoretical analysis on hydration thermodynamics of proteins J. Chem. Phys. 2006 125 024911 10.1063/1.2213980 40. Atkins P de Paula J Physical Chemistry 2006 8 Oxford Oxford University Press 41. Chiti F Dobson CM Protein misfolding, functional amyloid, and human disease Annu. Rev. Biochem. 2006 75 333 366 10.1146/annurev.biochem.75.101304.123901 16756495 42. Chiti F Stefani M Taddei N Ramponi G Dobson CM Rationalization of the effects of mutations on peptide and protein aggregation rates Nature 2003 424 805 808 10.1038/nature01891 12917692 43. Fernandez-Escamilla A-M Rousseau F Schymkowitz J Serrano L Prediction of sequence-dependent and mutational effects on the aggregation of peptides and proteins Nat. Biotechnol. 2004 22 1302 1306 10.1038/nbt1012 15361882 44. Tartaglia GG Vendruscolo M The zyggregator method for predicting protein aggregation propensities Chem. Soc. Rev. 2008 37 1395 1401 10.1039/b706784b 18568165 45. Maurer-Stroh S Exploring the sequence determinants of amyloid structure using position-specific scoring matrices Nat. Methods 2010 7 237 242 10.1038/nmeth.1432 20154676 46. Wang X Das TK Singh SK Kumar S Potential aggregation prone regions in biotherapeutics: A survey of commercial monoclonal antibodies mAbs 2009 1 254 267 10.4161/mabs.1.3.8035 20065649 47. Wang X Singh SK Kumar S Potential aggregation-prone regions in complementarity-determining regions of antibodies and their contribution towards antigen recognition: A computational analysis Pharm. Res. 2010 27 1512 1529 10.1007/s11095-010-0143-5 20422267 48. Buck PM Computational methods to predict therapeutic protein aggregation Methods Mol. Biol. 2012 899 425 451 10.1007/978-1-61779-921-1_26 22735968