Antibody association in solution: cluster distributions and mechanisms

ABSTRACT

Understanding factors that affect the clustering and association of antibodies molecules in solution is critical to their development as therapeutics. For 19 different monoclonal antibody (mAb) solutions, we measured the viscosities, the second virial coefficients, the Kirkwood-Buff integrals, and the cluster distributions of the antibody molecules as functions of protein concentration. Solutions were modeled using the statistical-physics Wertheim liquid-solution theory, representing antibodies as Y-shaped molecular structures of seven beads each. We found that high-viscosity solutions result from more antibody molecules per cluster. Multi-body properties such as viscosity are well predicted experimentally by the 2-body Kirkwood-Buff quantity, G₂₂, but not by the second virial coefficient, B₂₂, and well-predicted theoretically from the Wertheim protein–protein sticking energy. Weakly interacting antibodies are rate-limited by nucleation; strongly interacting ones by propagation. This approach gives a way to relate micro to macro properties of solutions of associating proteins.

KEYWORDS:

• Aggregation
• viscosity
• G22
• modeling
• thermodynamic perturbation theory

Introduction

Over the past several decades monoclonal antibodies (mAbs) have become one of the fastest-growing categories of therapeutics. mAbs are currently the most widely marketed biologic and the market value for mAbs is predicted to exceed $300 billion in 2025.¹ They are used to treat cancers, autoimmunity diseases, metabolic disorders and emerging infectious diseases. This broad range of applicability is due to their high binding affinity and specificity, long circulation half-life, and easy manufacturing.^1,2 Unlike most traditional drugs, biological therapeutics are delivered to patients as liquids.³ For mAbs to be effective for therapeutic use, they need to be administered to a patient in a relatively large dose. For subcutaneous injection, the maximum volume that can be delivered per dose in such cases is limited, thus requiring high mAb concentrations (e.g., as high as 175–200 mg/ml).⁴ High viscosity of antibody solutions at these concentrations causes difficulties in manufacturing and administration of these therapeutics.^4–6

In formulating biological drugs, it is desirable to create solutions that have high antibody concentrations and yet low viscosities, implying minimal association. We are interested in how the macroscale properties of liquid solutions of proteins, particularly monoclonal antibody (mAb) molecules, arise from their underlying intermolecular interactions. In several previous studies, it was observed that antibody association correlates with net hydrophobicity and hydrophobic patches,^7,8 with net charge on the proteins,^7–9 and with charge asymmetry.¹⁰ For example, by taking into account the hydrophobicity and charge of the amino acids in the protein sequence, Sankar et al. developed an algorithm called AggScore that identifies aggregation-prone regions in several well-studied proteins.¹¹ Mutational studies show that protein aggregation propensity correlates with positively charged surface residues (the more positive protein surfaces are, the less soluble the proteins are), the ratio of lysine to arginine content, and exposed hydrophobic patches.¹² While net charge matters, the distribution of charges, for example reflected in the dipole moment,¹³ particularly at high mAb concentrations is also important.^14,15

Past insights into the intermolecular interactions in protein associations have come from: (1) experiments on solution viscosities, liquid–liquid phase equilibria, second virial coefficients, and scattering structure factors^9,16–18; (2) atomistic molecular simulations (often by Molecular Dynamics, MD) using physics-based forcefields with appropriate solvation models^{16,17,19–22}; and (3) coarse-grained statistical mechanical solution theories.^23–27 Here, we combine coarse-grained Wertheim thermodynamic perturbation theory^28,29 with experimental data provided by Regeneron on 19 mAb systems. This Wertheim approach offers advantages over atomistic simulations, as it gives solution properties, such as viscosities, liquid phase properties and virial coefficients that are too computationally challenging for atomistic simulations. It also goes beyond experimental data alone in giving insights into cluster distributions and the energies and entropies that drive them. These quantities are otherwise difficult to determine. For example, for determination of cluster size distribution, experimental techniques, such as composition-gradient multi-angle static light scattering (CG-MALS),³⁰ and analytical ultracentrifugation³¹ are commonly used, but both require calibration standards, are limited with the detector resolution, the separation techniques are known to alter the nature of sample components, and the interpretation of the results is model dependent. Alternatively, attempts have been made to determine the cluster distribution in protein solutions using coarse-grained computer simulations.^{17,20,32–34} The computer simulations, however, are time-consuming, and their results depend on the protein model. Statistical mechanics on the other hand provides the basis for calculating the protein cluster size distribution quickly, and accurately. Below, we first describe our theoretical approach.

Methods

The 7-bead wertheim model of antibody association

Liquid-state statistical mechanics theories are recipes for computing macroscale solution properties from the microscale interactions of the underlying molecules. However, protein solutions pose two challenges for such approaches: (1) they often treat the solute (protein) as spheres, having minimal structure-property relationships and (2) such theories are usually low-density expansions dominated by few-particle interactions. Yet, proteins are often associated into multi-protein clusters. To address these two problems, we have recently^25,26 applied the Wertheim solution theory of strongly associating liquids²⁸ and we represent the molecular structure at a course-grained level; each antibody molecule has 7 beads forming a Y shape; see Kastelic et al.²⁷ (Figure 1). The model can then be treated with the associative Wertheim thermodynamic perturbation theory^28,29 to obtain the free energy, and related thermodynamic properties, as well as cluster size distribution.The latter can be, in combination with relations from polymer physics converted to the solution viscosity.^27,32 The model and theory have been successfully used to analyze the viscosity data of mAb solutions at different conditions, and in the presence of different excipients.^27,32 At the moment, however, the model, in-spite its successes, is missing its predictive value; the parameters for the model, as well as those for evaluating the contribution of different size aggregates to viscosity, were obtained by fitting the theoretical results to the experimental data. Details of 7-bead Wertheim theory are given in the SI.

Figure 1. 7-bead model of antibody associations.

Red and green patches are interaction sites through which molecules attract each-other via square-well attractive potential.²⁷ The gray patches represent infinite attractive interaction providing the hard spheres to be glued together, and forming a Y-shape molecule characteristic for antibodies. All the spheres are of the same size, and namely of diameter 2 nm to add up to the hydrodynamic radius of 5 nm as observed experimentally (R_H = 2.5σ). Following the experimental findings, molecules self-associate only through A and B interaction sites that are of the same strength (symmetric antibody molecule). u_KK represents the interactions strength between two k interaction sites on the protein surface that are separated by a distance r.

Here is how we apply the 7-bead Wertheim theory in the present work. First, we require an interaction distance (ω), which we take to be the same as in previous work (namely, ω = 0.18 nm which roughly corresponds to the length of a hydrogen bond.²⁷ Next, in full generality, the theory allows for three different types of interactions: AA, BB, and AB. Here, given that the data does not give sufficient granularity, we assume only a single type of energetic attractive interaction: ε_AA = ε_AB = ε₀. The attractive interactions between molecules being the cause for the self-association process leading to the formation of different size clusters in the solution, we have used the experimentally determined cluster-size distributions of the 19 mAb at low concentrations as a reference, and adjusted the interaction strength parameter to obtain the best agreement between theoretical and experimental data. This was shown previously to give excellent results for systems with strongly attractive directional forces, such as proteins and antibodies.^28,29

We are interested in the populations, P(n), i.e., the mass fractions of the different clusters containing n antibody molecules each. The results are shown in Figure S1, and the set of interaction parameters is given in Table S1. Then, from P(n) we compute the viscosities using¹²:

(1) $\mathit{ln}\mathit{\hspace{0.17em}}\left(\frac{\eta }{{\eta }_0}\right)=\mathit{\hspace{0.17em}}\sum _n^{\mathrm{\infty }}\mathit{\gamma }\mathit{P}\left(\mathit{n}\right)\mathit{c}{n}^d\mathit{\hspace{0.17em}}$(1)

where η₀ is the viscosity of the solvent, γ is the mass concentration of the solution, and where c, and d are constants describing how clusters contribute to viscosity. This model has been found to successfully describe the viscosities of antibody solutions of different pH, and with different salt concentrations, and different excipients.^17,20,27,32 In previous work different sets of c and d values have been used, obtained by fitting the calculated viscosity to experimental data. In the procedure, the interaction parameter, ε₀, and c and d values have been simultaneously modified, not considering the information regarding the experimental cluster distribution. In this work we have used the experimental cluster distributions from this work, as well as the ones obtained previously through combining SAXS measurements and computer simulations,^32,33 and determined the c and d parameters in the Equation (1)(1) $\mathit{ln}\mathit{\hspace{0.17em}}\left(\frac{\eta }{{\eta }_0}\right)=\mathit{\hspace{0.17em}}\sum _n^{\mathrm{\infty }}\mathit{\gamma }\mathit{P}\left(\mathit{n}\right)\mathit{c}{n}^d\mathit{\hspace{0.17em}}$(1) to properly transfer the cluster distributions to the experimentally determined viscosities. An important finding was that the c and d parameters are universal, not depending on the mAb, neither on the conditions in the solution. The values obtained were c=(0.017 ± 0.003)mL/mg, d=(0.5 ± 0.1), and allow a reliable conversion of cluster-size distribution to viscosities (± indicates the standard deviation throughout the text.). The values previously used to calculate the viscosity for a 7-bead model (c = 0.01205 mL/mg, and d = 0.3762)²⁷ also fall within the determined values.

To further interpret these values, we first rewrite Equation (1)(1) $\mathit{ln}\mathit{\hspace{0.17em}}\left(\frac{\eta }{{\eta }_0}\right)=\mathit{\hspace{0.17em}}\sum _n^{\mathrm{\infty }}\mathit{\gamma }\mathit{P}\left(\mathit{n}\right)\mathit{c}{n}^d\mathit{\hspace{0.17em}}$(1) in a form similar to Mooney’s well-known empirical viscosity equation for hard spheres³⁵ in which $\mathit{ln}\mathit{\hspace{0.17em}}\frac{\eta }{{\eta }_0}$ is proportional to ${\mathrm{\Sigma }}_n{\left[\eta \right]}_n{\mathrm{\Phi }}_n$, (${\left[\eta \right]}_n$ is the intrinsic viscosity of n-mer cluster, and Φ_n is its volume fraction):

(2) $\mathit{ln}\mathit{\hspace{0.17em}}\left(\frac{\eta }{{\eta }_0}\right)=\mathrm{c}\sum _{\mathit{n}=1}^{\mathrm{\infty }}{\mathit{\gamma }}_n{n}^d\mathit{\hspace{0.17em}}$(2)

where γ_n = γP (n) and represents the mass concentration of a cluster consisting of n monomers. Comparing the two expressions we can see that n^d plays the role of the intrinsic viscosity, [η], which is, for a n-mer cluster at low shear rates, written as³⁴ [η]_n ∝ n^{(3−df)/df.34} d_f is the so-called fractal dimension which has values between 1 (for linear clusters), and 3 (for compact spherical clusters),³⁶ corresponding to values of d between 0.5 (linear clusters), and 0 (spherical clusters). Our value 0.5 implies the presence of extended linear clusters in the solutions which is consistent with our model. Constant c, on the other hand, can be considered as a sort of proportionality constant, and is related to the intrinsic viscosity of the monomers. Our value 0.017 mL/mg is close to the experimentally obtained intrinsic viscosity, which has been found to be around 0.01 mL/mg for various mAbs.^37,38 The importance of the cluster shape for the viscosity of the mAb solutions has also recently been shown by Dandekar et al.³⁹

Results

This model captures well the experimental antibody viscosity curves

The solution viscosities computed by the 7-bead Wertheim theory with the experimental data on Regeneron’s 19 antibody systems are shown in Figure 2.

Figure 2. Viscosities as a function of antibody concentrations. In all cases the lines represent the results for the model, and the symbols were obtained experimentally. For each sample, the viscosity has been measured 10 times, the results given are arithmetic averages, the error bars approximately corresponding to the size of the symbols. All results applied to mAbs in pure water at 293 K.

For all 19 mAbs, an excellent agreement between experimental and theoretical viscosities is obtained, which enabled us to use the information about the cluster distribution from the theory for further insights into the mechanism of cluster formation in mAb solutions.

Bigger clusters contribute more to the viscosity

Our modeling predicts that different antibodies have different distributions of cluster sizes (Figure 3; see the Supporting Information (SI) for all 19 systems).

Figure 3. Cluster size distribution vs. protein concentration. Theoretical predictions for normalized mass fraction distribution of clusters of size n, P(n) for our model as a function of concentration for two mAbs (A – left, and H – right). The results for all mAbs studied here are given in SI (fig SI 2).

Two antibodies, A and H, show the paradigmatic extremes among the 19 systems we studied. Although mAbs A and H belong to the same antibody isotype, IgG1, A-like molecules form mostly monomers and dimers over the full concentration range, whereas H-like molecules form higher-order clusters at higher concentrations. Both scenarios occur regardless to the mAb class, IgG1 (mAbs A, F, and H, and RP7, RP8, RP9, and RP10), or IgG4 (mAbs B, C, D, E, G, I, and J, and RP1, RP2, RP3, RP4, and RP6), confirming predictions that the Fc region of mAbs plays a minor role in protein–protein interaction. Confirmation of this observation lies in the fact that Fc regions of the remaining mAbs do not differ significantly from each other, but nevertheless these antibodies clearly show different viscosities.

Figure 4 shows the more granular breakdown from the model of how much each of the cluster types contributes to the viscosity (all 19 are given in Figure S3). All results of cluster sizes vs. viscosity simply reflect the intermolecular interaction strength (weaker in A-like antibodies, stronger in H-like antibodies).

Figure 4. The relative contributions of different-sized clusters to viscosity. Predictions of the wertheim model for the relative contributions of protein monomers, dimers, 5-mers and 10-mers to the solution viscosities, vs protein concentration, for mAb A, and mAb H. The contributions to the viscosity of other size clusters (trimers to 9-mers, and higher order oligomers) are not shown in the figure.

Cluster-size distributions reflect their kinetic assembly mechanism

Roberts et al. have shown how cluster-size distributions in protein solutions can be related to their kinetic mechanisms of cluster formation through the Lumry – Eyring nucleated polymerization theory.^40,41 In short, the rate limiting step to cluster formation can either be the nucleation step or the elongation step. According to this theory, the aggregate growth mechanism can be determined from a plot of the weight average molecular weight, including monomers and aggregates in solution, M^tot/M^mon vs (1−m)² (m is the fraction of monomers, compared to the initial value, M^tot is the total weight-average molecular weight, w and M^mon is the monomer molecular weight.⁴¹ When this plot is linear, it means aggregates grow primarily through monomer addition to an already existing aggregate (the so-called chain polymerization (CP) case). Otherwise, the mechanism is nucleation dominated (ND) and the aggregates do not reach sizes larger than dimers and small oligomers. An upturn in these types of plots represents the growth of aggregates by aggregate – aggregate condensation polymerization (AP).⁴¹ Figure 5 shows Lumry – Eyring plots, computed in the 7-bead Wertheim model, for antibodies A and H (remainder are shown in Figure S4).

Figure 5. Lumry – Eyring plots for learning the rate-limiting steps of formation, for mAb a (left) and mAb H (right) as obtained by our model treated with wertheim perturbation theory. A forms by the ND mechanism, while H forms from the CP mechanism.

For mAb A, the clusters form directly from monomers (ND), so dimers prevail in the solution. For mAb H, the clusters form by elongation (CP), where new monomers preferentially bind to already formed chains rather than to other monomers. By analyzing the data for all 19 mAbs studied in this work we found that the nucleation dominant mechanism gives lower viscosities, while elongation mechanism gives larger viscosity increases with protein concentration (Figure 6). One implication is that viscosity could be reduced by energetically or sterically restricting one Fab arm, which would limit its ability to propagate in growing chains.

Figure 6. Solutions of lower viscosities are nucleation-rate-limited; higher viscosities are propagation limited, as determined from Lumry – Eyring plots. For mAbs A, B, F, C, RP3, RP7, and RP8, nucleation is the dominant mechanism for cluster formation. For mAbs G, D, I, J, H, E, RP1, RP2, RP4, RP6, RP9, and RP10, elongation following initial dimerization is rate limiting.

Discussion

Predicting concentrated solutions from dilute-solution measurements

The ability to interpret multi-body properties of colloid solutions in terms of 2-body properties that can be measured or modeled more easily and in dilute solutions has long been desired. Could we predict liquid-phase equilibria and viscosities, for example, from simply knowing pairwise protein interactions in dilute solutions through measuring the second virial coefficient, B₂₂, or obtained from the scattering of light or neutrons, or from equilibrium analytical ultracentrifugation, self-interaction chromatography, or osmotic pressure experiments.^22,42–46B₂₂ can be expressed as an integral over solute–solute pair distribution function, g₂₂, in the limit of low solute concentration⁴⁷:

(3) $B_{22}=\mathit{\hspace{0.17em}}{B}_{22}^{HS}-\mathit{\hspace{0.17em}}\frac{1}{2}\int \left(g_{22}\left(\mathit{r}\right)\mathit{\hspace{0.17em}}-\mathit{\hspace{0.17em}}1\right)\mathit{\hspace{0.17em}}4\mathit{\pi }{r}^2\mathit{\hspace{0.17em}}\mathit{dr}\mathit{\hspace{0.17em}}=\mathit{\hspace{0.17em}}{B}_{22}^{HS}-\mathit{\hspace{0.17em}}\frac{1}{2}\int \left(e^{-{\mathit{u}}_{22}/k_B\mathit{T}}-\mathit{\hspace{0.17em}}1\right)\mathit{\hspace{0.17em}}4\mathit{\pi }{r}^2\mathit{\hspace{0.17em}}\mathit{dr}$(3)

where u₂₂ is the orientation-averaged pair potential between two solute (protein) molecules. Negative B₂₂ values indicate a net attraction while positive values indicate protein-protein repulsion.⁴³ However, B₂₂ is sometimes considered a surrogate simple measurement that could predict higher-concentration behaviors, a more principled quantity for higher concentrations is the corresponding Kirkwood – Buff integral, G₂₂:

(4) $G_{22}=\int \left(<g_{22}\left(\mathit{r}\right)>-\mathit{\hspace{0.17em}}1\right)\mathit{\hspace{0.17em}}4\mathit{\pi }{r}^2\mathit{\hspace{0.17em}}\mathit{dr}=\int \left(e^{-{\mathit{w}}_{22}/k_B\mathit{T}}-\mathit{\hspace{0.17em}}1\right)\mathit{\hspace{0.17em}}4\mathit{\pi }{r}^2\mathit{\hspace{0.17em}}\mathit{dr}$(4)

where u₂₂ is replaced by the potential of mean force, W₂₂, and <g₂₂(r)> is the average molecular pair correlation function.^{9,13,19,42,43,48,49} Similarly as B₂₂, negative values of -G₂₂ indicate a net attraction while positive values indicate net protein–protein repulsion. As such, the protein–protein Kirkwood-Buff integral is a measure of the net protein–protein interaction strength in the solvent.⁵⁰ In short, the difference is that B₂₂ describes two proteins interacting in the absence of solvent while G₂₂ describes the interactions in the presence of solvent.

Figure 7 (and Table S2) tests whether either of these two pairwise quantities, B₂₂ or G₂₂, is a predictor of solution viscosities for our set of 19 antibodies at the mAb concentration of 175 mg/mL. The results shown here are clear: G₂₂ is an excellent predictor of concentrated solution viscosities, while B₂₂ is not correlated at all. This result is consistent with earlier work ^{9,13,19,42,48,49} and indicates the importance of using the Kirkwood-Buff integral, G22, which can be determined by static light scattering¹³ or SAXS/SANS measurements.³³ This is consistent with the findings of Ghosh et al.,⁴⁹ and Barnett et al.,⁵¹ who found G22 may offer semiquantitative means to predict aggregation mechanism

Figure 7. Left: B₂₂ is not predictive of solution viscosities. Right: G₂₂ is predictive. Correlations between experimental viscosities at mAb concentration 175 mg/mL. The colors indicate ranges of B₂₂ that are experimentally considered for classification of colloidal stability of antibodies: green – repulsive, yellow – near ideal, red – attractive. For each sample, B₂₂ has been measured in triplicates, the error bars approximately corresponding to the size of the symbols.

To test the possible importance of G₂₂, as a predictor of solution viscosities, the coefficient was determined experimentally for all 19 mAbs under investigation at different mAb concentrations. We determined that the G₂₂ remains approximately constant above mAb concentration 80 mg/mL (the results are for concentration 80 mg/mL given in Table S2, and the concentration dependence in Figure S5). This indicates that the net protein–protein interaction remains constant beyond this concentration. From Figure 8 it can be seen that, as predicted by our theory, the viscosities at high mAb concentration are correlated to G₂₂, indicating the importance of the mAb-mAb interactions that are beyond pair-wise ones.

Figure 8. The correlation between experimental viscosities at mAb concentration 175 mg/mL and experimentally determined G₂₂ coefficients at mAb concentration 80 mg/mL. For each sample, G₂₂ has been measured in triplicates, the error bars approximately corresponding to the size of the symbols. The linear fit (red line) is represented by the equation y = −92.25x + 21.66; Pearson’s r = −0.65.

What types of interaction dominate antibody associations in solution?

Our 7-bead Wertheim model has one parameter, ε₀, representing the protein–protein short-ranged attraction energy between the particular sites indicated in Figure 1. First we note, as Figure 9 shows, that predicted solution viscosities correlate quite well with this interaction strength. In particular, it can be well fit using the relationship: ln(η(175 mg/mL)/cP) = −4.902 + 0.9588 (ε₀/(kcal/mol)) (similar relationships can be obtained for other mAb concentrations).

Figure 9. Solution viscosities correlate with, ε₀, the protein-protein affinity of the model at mAb concentration 175 mg/mL. Pearson’s r for the right panel is 0.99.

In this way, the Wertheim theory reduces the prediction of macroscopic solution behavior to the microscopics of the intermolecular interactions through a single fitting parameter, the attraction strength. Interestingly, the above relation resembles the shear viscosity model where viscosity is a thermally activated process in which a molecule, in order to move to a neighboring free space, has to overcome activation energy barrier, E_a, created by the resistance of the surrounding building units (ln(η) = lnA_s + E_a/RT).^42,52,53 Within this theory our association parameter reflects the microstructure of the solution, its configurational entropy in particular (e.g. the configurational entropy of the solution is inversely proportional to activation energy,⁵⁴ and could in principle be determined from the temperature dependence of the solution viscosity. We will investigate this further in our future work.

Even though the interaction parameter of our model only describes an average local interaction, and as a single parameter cannot comprise different types of interactions known to define the colloidal stability of protein solutions, the parameter was found to be directly correlated to the high concentration solution viscosity (Figure 9). Several studies exist where this quantity has been predicted from the above mentioned mAb properties. We therefore examined the possible correlations between the physical properties of the 19 mAbs studied here, and the interaction parameter of our model; details are given in SI. Even though previous studies that combine homology modeling with different machine learning techniques have found correlations with mAb molecule properties, such as hydrophobicity, charge symmetry, and net charge,^{2,9,10,12,55–58} no conclusive correlations are observed here, even when each mAb isotype is evaluated separately (see SI for details).

Conclusions

We studied the association of monoclonal antibody molecules in liquid solutions, through comprehensive experiments on 19 systems combined with a 7-bead liquid statistical mechanics theory. Some antibodies are stickier than others. The stickiest ones tend to associate into larger clusters, leading to disproportionately higher-viscosity solutions. Viscosities and other multi-body properties can be anticipated quantitatively through experiments on simpler protein pairwise properties in dilute solution, such as the Kirkwood-Buff property, G₂₂. Correspondingly, the sticking energy quantity in the theory is also directly, but nonlinearly, predictive of multi-body properties. Finally, when combined with the Lumry-Eyring theory, we can predict the kinetic cluster formation mechanism. Weak binders tend to be nucleation-rate-limited and strong binders are propagation-rate-limited.

Supporting information

Supporting information consists of details on experimental procedures, additional theoretical background, additional figures of our detailed mAb analysis and also results of homology modeling.

Acknowledgments

S.B., K.A.D., and B.H.-L. acknowledge the support from the National Institutes of Health (NIH) RM1 award “Solvation modeling for next-gen biomolecule simulations” (grant No. RM1GM135136).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/19420862.2024.2339582

Additional information

Funding

The work was supported by the National Institutes of Health [RM1GM135136].