Quantification of natural abundance NMR data differentiates the solution behavior of monoclonal antibodies and their fragments

ABSTRACT

Biotherapeutics are an important class of molecules for the treatment of a wide range of diseases. They include low molecular weight peptides, highly engineered protein scaffolds and monoclonal antibodies. During their discovery and development, assessments of the biophysical attributes is critical to understanding the solution behavior of therapeutic proteins and for de-risking liabilities. Thus, methods that can quantify, characterize, and provide a basis to inform risks and drive the selection of more optimal antibody and alternative scaffolds are needed. Nuclear Magnetic Resonance (NMR) spectroscopy is a technique that provides a means to probe antibody and antibody-like molecules in solution, at atomic resolution, under any formulated conditions. Here, all samples were profiled at natural abundance requiring no isotope enrichment. We present a numerical approach that quantitates two-dimensional methyl spectra. The approach was tested with a reference dataset that contained different types of antibody and antibody-like molecules. This dataset was processed through a procedure we call a Random Sampling of NMR Peaks for Covariance Analysis. This analysis revealed that the first two components were well correlated with the hydrodynamic radius of the molecules included in the reference set. Higher-order principal components were also linked to dynamic features between different tethered antibody-like molecules and contributed to decisions around candidate selection. The reference set provides a basis to characterize molecules with unknown solution behavior and is sensitive to the behavior of a molecule formulated under different conditions. The approach is independent of protein design, scaffold, formulation and provides a facile method to quantify solution behavior.

KEYWORDS:

• Nuclear magnetic resonance spectroscopy
• biophysical characterization
• antibody characterization
• antibody and antibody-like molecules
• singular value decomposition

Introduction

Protein-based therapeutics like monoclonal antibodies (mAbs) have found important clinical applications for a myriad of diseases. Their approval for clinical use has increased dramatically, with over 30 new mAbs being approved since 2018.^1–3 The traditional IgG mAb scaffold adopts a canonical fold in which two heavy and two light chains associate to create two antigen-binding fragments (Fab), composed of the variable heavy (VH) and variable light chain (VL) that can bind to a target.⁴ The Fabs are attached by a flexible hinge to the Fc domain, which is composed of heavy chain constant regions (CH2, CH3) that dimerize. In addition, the conventional mAb scaffold can also be modified by mutation to elicit unique functional behavior, including increasing half-life,⁵ organized higher-order multimer structure formation⁶ and reduced Fc receptor engagement,⁷ but other modifications have been documented.^8–11 The clinical success of mAb therapeutics has spurred innovation and the engineering of a multitude of different therapeutic proteins that engage targets using unique structural folds. Therapeutic proteins cover a spectrum of scaffolds that range from small single variable heavy chain domains (approximate molecular weight (MW): 12–15 kDa), typically referred to as single-domain V_HH monomers,¹² to self-associating mAbs (MW: 100s – 1000 kDA hexameric IgG, IgM, respectively).^6,13 Some of these scaffolds or their fragments, can also be tethered to form multidomain proteins that are referred to as multivalent molecules when they bind two or more antigens.¹⁴ Tethering allows creation of molecules that can engage multiple different epitopes on the same target or cross-link different targets. Furthermore, antibody-drug conjugates also present an exciting method to deliver payloads to specific cells.^15,16 The availability of different protein scaffolds presents an exciting opportunity to discover novel therapeutics with differentiated mechanisms of action.^17,18 However, antibody and antibody-derived protein therapeutics can present a developability challenge that can only be overcome by leveraging biophysical and structural methods that characterize, quantify and differentiate their solution behavior.^19,20 Novel methods that can help differentiate between candidates using simple comparators that are sensitive across different formulation conditions are thus needed.

The biophysical characterization of therapeutic proteins serves a critical role in guiding decisions within discovery and developability efforts.^20–24 Chromatography, electrophoresis, thermal-melt fluorimetry, mass spectrometry and light-scattering techniques largely anchor many platforms that are used to qualify proteins that have therapeutic potential.^25–29 These approaches can capture the average state of the protein, informing on size, charge and distribution of aggregated states. Some of these approaches are also sufficiently sensitive to detect site-specific modifications that can manifest as liabilities.²⁶ Generally, these techniques are combined to develop fit-for-purpose platforms to qualify protein behavior. Over the past few years, Nuclear Magnetic Resonance (NMR) spectroscopy has emerged as a useful tool to profile a range of mAbs and other biologics. NMR is a technique that can profile molecules in solution, independent of formulation, while maintaining atomic resolution.^30,31 NMR can also be sensitive to species that are sparsely populated (>1%) as long as they maintain favorable NMR relaxation properties.³² One-dimensional (1D) NMR-based methods are information rich and have been used to profile mAb and mAb-like molecules.³³ Herein, we will denote the collection of different classes of mAb, mAb-like and mAb fragment molecules, referred to here as Ab(L/F)s. The 1D spectrum can report on the overall fold and, with careful normalization, molecules with a common overall structural topology can be compared.³¹ In unique situations, if particular resonance peaks can be identified, site-specific modifications can be monitored for changes over time and detect different effects that affect the molecule studied.³¹ However, spectral overlap can still be a limitation with detailed analyses.

The application of Diffusion-NMR (commonly referred to as Diffusion Ordered Spectroscopy (DOSY) NMR) has also been demonstrated to be useful in directly measuring the translational diffusion rate of biological molecules.³⁴ An advantage of DOSY NMR is the ability to profile molecules that are <1 kDa to several MDa, and it is sensitive to the distribution of translational diffusion rates that could be sampled by any detectable species.³⁴ This can enable the quantification of stably formed aggregates. Coupled with¹H NMR relaxation-based approaches, this provides additional ways to probe the average differences within mAb and mAb-like molecules.^34,35 Modern NMR instrumentation has dramatically improved the sensitivity of conventional NMR experiments. This breakthrough, coupled with sparse sampling techniques, has enabled the detection of two-dimensional (2D) [¹H, ¹³C]-HSQC or [¹H, ¹³C]-HMQC fingerprint spectrum (also referred to as a methyl fingerprint spectrum), which report on the behavior of individual terminal methyl and/or aliphatic groups within mAbs and Fabs at natural abundance.³⁰ This obviates the necessity to isotope label the samples using recombinant or cell-free expression systems. Although, examples of Ab(L/F)s being isotope labeled have been documented, this process is very challenging and not accessible to the many targets and scale typically required in a discovery setting.^31,36,37 Natural abundance 2D NMR applied to mAbs and Fabs has largely been championed by the Marino group, where they demonstrated multiple applications including detected differences in mAb glycosylation, reported a statistical approach to assess large NMR signal differences between a set of biologics samples, and established the robustness and variability of 2D [¹H, ¹³C]-HSQC fingerprinting for biosimilar comparison.^{30,36,38–41} However, extracting and comparing quantitative aspects between different types of Ab(L/F)s still remains a challenge.

A collection of 2D [¹H, ¹³C] correlation spectra for a variety of Ab(L/F)s, with orthogonally qualified solution behavior, would provide a basis of comparison against novel molecules. We propose that the linewidths and relative peak intensities of the resonances from 2D correlation spectra should be used as a measure of contrasting solution behavior between different Ab(L/F)s. This is because the NMR linewidth is a direct reporter on the overall rotational tumbling time of the molecule and conformational exchange events.³² The overall rotational tumbling time is directly related to the shape/size of the molecule and is the primary source contributing to relaxation phenomena, which occurs on the nanoseconds timescale. Conformational exchange manifests as an addendum to the NMR linewidth and reports on lowly populated transient processes that occur on the micro- to millisecond timescale. These processes can involve self-association and intra-/interdomain fluctuations that can induce unique structural states.^42–44

From the above summary, linewidths and relative peak intensities, from different Ab(L/F)s could be used as proxies for NMR relaxation rates. Here, we propose that large random sampling of a subset of these spectral parameters and their incorporation into a Principal Covariance Analysis (PCA) simplifies their interpretation. This analysis provided a simplified basis to compare different Ab(L/F)s, and a trend emerged in which the largest components of the input data directly correlate with molecular hydrodynamic properties. By using a collection of different Ab(L/F)s, we could test the behavior of molecules with unknown behavior. The analysis also shows a sensitivity to changes in construct design that contribute to decisions around candidate selection. Furthermore, it is shown that the approach is sensitive enough to readout changes associated with a mAb formulated under three different conditions.

Results

Random sampling of principal components reveals an excellent correlation with hydrodynamic properties of Ab(L/F)s

We acquired natural abundance 2D [¹H, ¹³C]-HMQC datasets for 20 different molecules that span four different Ab(L/F) classes with MWs between 12 and 150 kDa (Table S1). Figure 1 displays examples of 2D NMR datasets from the four different classes of Ab(L/F)s evaluated here. The set included individual V_HH domains (Figure 1a; red), multivalent or tethered V_HH domains (Figure 1b; blue), antigen-binding fragment (Fab; Figure 1c; green) and full-length mAb (Figure 1d; cyan). From the 2D NMR datasets, only qualitative differences between the classes of molecules are evident in Figure 1a-d. The classes were specifically chosen to enable construction of a set that could distinguish features based on their modular construction. As a rough approximation, the selected classes are composed of single to multiples of small immunoglobulin domain(s) that are folded or tethered together.

Figure 1. (a-d) [¹H, ¹³C]-HMQC correlation spectra for four different classes of Ab(L/F)s recorded at natural abundance. The four different classes of Ab(L/F)s that are represented include an individual V_HH domain (a; red), multivalent V_HH domains (b; blue), a fragment antigen binding domain (Fab) (c; green) and a full monoclonal antibody (mAb) (d; cyan). The [¹H, ¹³C]-HMQC datasets reveal resonance peaks which correspond to individual methyl groups and reveal that high quality structural fingerprint spectra of Ab(L/F)s (with molecular weights (MW) spanning approximately 12–150 kDa) can be collected without the need of isotope labeling. As we progress from the smaller MW Ab(L/F)s (a, b) to the large fAb (c) and mAb (d) particles the number of peaks can vary and the “relative sharpness” of the peaks change significantly. The change in peak quality is related to its intensity and linewidth, which can be quantitated by the apparent linewidth (LW_app) for both the ¹H and ¹³C dimensions. (e) No correlation exists between LW_app for both ¹H and ¹³C dimensions across individual peak resonances and is largely indistinguishable between different types of protein therapeutics. The plot has the LW_app plotted for twenty different protein therapeutic molecules. The data points are colored with respect to the type of protein therapeutic as defined in panels a-d. (f) Low correlation is observed between the average LW_app between ¹³C and ¹H for twenty different Ab(L/F)s. Each data point represents the average LW_app and the error bars are derived from the standard deviation across all values for a given Ab(L/F) sample

These datasets report on individual methyl-bearing amino acids and some methylene groups within each molecule. The smaller molecules (Figure 1a, b) maintain a smaller number of peaks as compared to larger counter parts because Fab and mAb molecules (Figure 1c, d) contain more methyl groups. An important difference between all class of molecules is the heterogeneity in peak intensities and resonance peak linewidths in both ¹H and ¹³C dimensions. For each observed resonance, within a 2D spectrum, two linewidths can be extracted corresponding to the ¹H nuclei and its bonded ¹³C atom. The linewidth of a resonance peak is directly associated to NMR relaxation phenomenon that are modulated by different physical processes that affect the molecule.³² Many approaches have been developed to accurately measure relaxation rates, which can be analyzed to determine parameters that define molecular reorientation and conformational exchange.^32,44 However, with Ab(L/F)s at natural abundance, the accurate extraction of these rates in an artifact-free fashion would require unrealistically long acquisition times for small molecules. For a majority of larger Ab(L/F)s, this would also not be amenable because of signal-to-noise restrictions. It is important to note that linewidths from [¹H, ¹³C]-HMQC datasets report on multiple relaxation pathways because of the generation of anti-phase and in-phase coherences during chemical shift labeling in both direct (¹H) and indirect (¹³C) dimensions, and are not directly comparable to NMR-based methods, which extract precise relaxation rates.⁴⁵ Therefore, we refer to linewidths extracted from a [¹H, ¹³C]-HMQC spectrum as an apparent linewidth (LW_app). This complexity can also complicate procedures that analyze line-shapes of resonance peaks, unless time-dependent changes in coherence evolution are taken into account.⁴⁶ However, as long as experimental acquisition parameters are constant, across all datasets, comparison is possible. Therefore, we propose that linewidths and relative intensities of resonance peaks detected from [¹H, ¹³C]-HMQC (Figure 1a-d) can serve as a proxy to report on the aforementioned physical processes and can provide a basis to quantitatively compare different Ab(L/F)s.

An automated procedure is required to enable facile analysis of natural abundance [¹H, ¹³C]-HMQC spectra of antibody and antibody-like molecules. In order to remove human bias from the peak selection protocol, [¹H, ¹³C]-HMQC spectra we used the automatic peak selection tool from the NMRPipe software package,⁴⁷ where resonances were selected based on a minimum signal-to-noise ratio criteria and were identified from within a constant spectral window (see Materials and Methods). A comparison of the linewidths from detected resonance peaks rendered no clear pattern that would enable quantitative comparison amongst different antibody and antibody-like molecules (Figure 1e, f). The raw individual LW_app values make it difficult to draw meaningful conclusions about molecular behavior. However, we note that the single V_HH domain resonances maintain smaller ¹³C LW_app as compared to the larger molecules (Fab and mAb, green and cyan Figure 1e, 1f), and when the average LW_app, in both dimensions, is compared (figure 1f), a qualitative binning is achievable. In general, the average LW_app for smaller Ab(L/F)s, like V_HH domains, show some separation between larger molecules. However, within a specific class the dispersion among the data points is poor and highly overlapped. Therefore, quantitative comparison by bulk LW_app parameters is of very limited value. In order to provide a basis, for which 2D spectra of Ab(L/F)s at natural abundance can be directly compared, we developed and evaluated a numerical approach that could highlight changes within spectral parameters.

An approach that can uniformly “normalize” the data from any 2D spectrum could allow the detection of unique differences, potentially providing the ability for more quantitative comparison between Ab(L/F)s. We hypothesize that using peak parameters from an array of different molecules, which maintain different structural scaffolds, could allow the identification of subtle differences between molecules. A variety of numerical approaches have been developed to assign quantitative similarities or dissimilarities between datasets that are complex in nature or that are dependent on many variables.⁴⁸ In particular, Principal Component Analysis (PCA) is a powerful approach that can take the same type of data from any number of samples composed of any number of observables that are dependent on any number of variables and can reduce their complexity by projecting their components into a simplified representation that allows facile detection of correlated or uncorrelated aspects between the input datasets.⁴⁸ Data subjected to a PCA are simplified, whereby the number of random variables that contribute to the data are reduced by projection into a lower dimensionality space. This projection identifies the minimum number of components that contribute to the input datasets and can establish how they might be distinguishable from other ones used in the calculation (i.e., how much is the variance between them). The power to distinguish subtle differences that are not readily apparent from the original data has led to PCA being used in a myriad of applications.^49,50 PCA approaches have also been applied to 2D NMR spectra in cases when the atomic peak assignments are accessible. For example, PCA has been applied to ascertain complex-binding mechanisms during NMR-based titrations where changes to one assigned polypeptide sequence are monitored upon step-wise addition of a ligand.⁵¹ This technique can also be used to evaluate the behavior of one molecule that is subjected to different solution or instrument conditions.^38,40 One criterion for PCA is that the input data array, across all samples, must have an equivalent number of values. Therefore, the direct application of PCA to 2D NMR, of different Ab(L/F)s, is not straightforward because the structural scaffolds, size and behavior of samples will generate a different number of resonance peaks within each spectrum. Here, peak assignment of many Ab(L/F) samples is not readily feasible, as uniform isotope labeling is not always possible, and/or the MW of the molecule can be large making assignment experiments intractable. Furthermore, we aim for any approach to remain indifferent to the type of molecules that could be tested, and thus, any method would benefit if it could operate without the need to label.

In order to avoid the necessity of peak assignment and enable comparisons of the natural abundance [¹H, ¹³C]-HMQC spectrum of any Ab(L/F), we combined the PCA with a procedure in which a subset of resonance peaks is selected, at random, from each spectrum and their PCs are calculated (see Materials and Methods). This random sampling of spectral information is repeated in a Monte Carlo fashion and the distributions of the calculated eigenvectors that arise from the repeated PCA are analyzed. We have termed the procedure “A Random Sampling of NMR Peaks for Covariance Analysis”. Hereafter, we will refer to this procedure as a Covariance Analysis (CA) method. From each spectrum, the LW_app in both ¹H and ¹³C dimensions and normalized peak amplitude (normalized to maximum peak intensity) were used as input data for this procedure.

All molecules were processed through this CA and the top PCs were analyzed (Figure 2). The first four PCs account for approximately 90% of the entire reference dataset and remain constant with limited variability throughout different iterations of the calculation (Figure S2). The values of the first two PCs that were determined from this CA method showed a separation between the different classes of antibody and antibody-like molecules. We hypothesized that since the input data are proxies for NMR relaxation phenomena, which is largely affected by the size and shape of the molecule,³² the first two PCs could correlate with parameters that are sensitive to this attribute. To that end, we determined the apparent hydrodynamic radius (R_h) for all the reference molecules using Diffusion NMR in order to compare with the results from the CA³⁴ (see Materials and Methods for details regarding the Diffusion NMR data). We compared the magnitude of the eigenvectors that were extracted for the reference set of antibody and antibody-like molecules and their MWs (Figure 2a, 2c) and to their R_h (Figure 2b, 2d). Interestingly, a very high correlation is observed between the first two PCs and these macroscopic properties with Pearson correlation coefficients (r) between 0.81 and 0.89. This demonstrates that peak parameters that are extracted from 2D NMR datasets and that are processed through this CA method are sensitive to properties of molecules. However, MW does not completely capture the distribution of values that came from this method. Although, the Pearson correlation coefficients of PC1 and PC2 between MW and R_h, are nearly equivalent, a reduced variance is observed when PC1 and PC2 are compared to the apparent R_h of the molecule (Figure 2b and 2d). This is because R_h is sensitive to not only the size, but the overall shape of the molecule, which is dependent on the tertiary and quaternary structure of the molecule. This approach also allows us to generalize features from the reference dataset. The red dashed lines in Figure 2 are linear fits to the data and permit back calculation of an observable parameter like R_h. The variance for larger molecules like multidomain V_HH (Figure 2; blue points), Fab (Figure 2; green points) and mAb (Figure 2; cyan points) molecules are reduced. However, the variance remains high in both PC1 and PC2 for individual V_HH domains when they are compared to MW and R_h. This indicates that there may be other effects that contribute to PC1 and PC2 for individual V_HH domains that are not captured by R_h. Our numerical approach is an application of PCA that can quantify natural abundance NMR data from a range of different Ab(L/F)s and reveals how relative shape and size can be discerned and rank ordered.

Figure 2. The molecular weight (MW) (a, c) and hydrodynamic radius (R_h) (b, d) of an array of Ab(L/F)s correlate very well with the first Principal Components (PC) derived from natural abundance 2D correlation spectra. All data points within the plots were generated from the CA protocol involving twenty different natural abundance 2D datasets. The points are colored with respect to the class of Ab(L/F). They include V_HH domains, multivalent V_HH domains, fragment antigen binding (Fab) domains, and monoclonal antibodies and are colored red, blue, green and cyan, respectively. A comparison of the MW and R_h with respect to the first two Principal Components (PC) reveals a strong correlation. Although a high Pearson correlation coefficient (r) is observed for MW it is improved when R_h is used as it is a better measure of the shape/size of a molecule in solution

Visualizing the first two PCs and the effect of higher order principal components

Our CA protocol can be applied to Ab(L/F)s with unknown solution behavior. By having a well profiled reference set, a test molecule, with unknown behavior can be benchmarked against previously characterized molecules (vide infra). However, within the reference set of molecules tested here, a comparison of PC1 and PC2 already highlights that some molecules within the same class can have unique features (Figure 3a). Directly plotting PC1 against PC2 (Figure 3a) constructs a “molecular ruler”; a measurement that is very sensitive to subtle differences due to the inclusion of multiple types of molecules in the analysis. The plot is linear in nature, indicating correlation between the largest contributor to LW_app, the molecule’s overall tumbling time, and R_h. The linearity is due to the fact that the MW of proteins, between 10s to 100s kDa, will cause LW_app to scale $L{W}_{app}\propto {\tau }_c/\left(1+{\left(\omega {\tau }_c\right)}^2\right)$; τ_c and ω are the overall tumbling time and the Larmor frequency for a given type of nucleus. This presentation provides a facile way to evaluate any additional molecule or to identify a unique feature within a molecule, when they are compared to one another. For example, one of the mAbs included in the reference set was the NISTmAb,⁵² a well-studied biologics reference standard. Interestingly, the NISTmAb (cyan point highlighted by the black arrow labeled by an asterisk (Figure 3a)) falls within a region where smaller Ab(L/F)s are located and is distinct from the other mAb molecules. This result agrees with reports of NISTmAb R_h measurement (by Dynamic Light Scattering) of 4.9 nm, which is lower than the 5.1 to 5.5 nm range for typical IgG mAbs.^53,54 This is distinct from the Fabs that were included in the reference set (Figure 3; green points). The Fabs largely cluster together indicating that, within these select molecules, they share a similar R_h. A similar trend is observed for the individual V_HH domains (red points; Figure 3a), but with improved separation between them as compared to direct evaluation against MW or R_h (Figure 2).

Figure 3. PC3 and PC4 do not report the same correlation as compared to PC1 and PC2. PC3 and PC4 are also sensitive to the impact of interdomain interactions between multivalent V_HH molecules. (a) Three multivalent V_HH molecules (highlighted with arrows in the second subset) consist of the same individual V_HH domains but with different linker lengths. The data point highlighted with an arrow attached to an asterisk corresponds to the NISTmAb (b) A comparison of PC3 and PC4 show a deterioration away from the correlation with R_h (Figure 1). Arrows indicate two multivalent V_HH molecules whose linker lengths differ. (c) The bottom subset represents a rank-ordering of the magnitude from PC3. A general trend appears where certain classes of Ab(L/F)s can be separated by the magnitude of PC3. This may indicate that in general Fabs, multivalent V_HH domains, individual V_HH domains, and mAbs larger to smaller inter-domain fluctuations, respectively. However, this approach is only sensitive to dynamic fluctuations that can manifest on a timescale that can significantly impact a resonance’s linewidths. The absence of significant value to PC3 for mAbs could indicate that their fluctuations may occur on another timescale (faster than μs or longer than ms).The blue point with an asterisk is a multivalent V_HH domain that has behavior similar to single V_HH domains. (d-f) The sensitivity of this CA to detect differences of Test-mAb formulated under different conditions. In d-f the natural abundance spectrum of Test-mAb within a different condition (Test-mAb d-f) was added as an additional dataset and compared against the twenty original Ab(L/F) datasets. The average apparent behavior of Test-mAb formulated under different conditions can be distinguished based on a comparison of where the dataset lies with respect to the plot of PC1 and PC2

When PC1/PC2 are compared, increased scatter is observed for the multivalent V_HH molecules (blue points; Figure 3a). This class of Ab(L/F) can be tethered by engineered amino acid linkers, and the effect of linker composition and length with respect to solution conformation of the whole molecule remains poorly understood.^15,18,55 It is important to have tools that can isolate differences in a molecule’s valency and linker length. This can affect the presentation of the antigen-binding surface, which can affect affinity, but also introduce aberrant intermolecular interaction that could lead to aggregation or stability issues.^56–58 In particular, understanding the effect of linker length can be empirical and relies on orthogonal in vitro datasets that could provide degenerate results. We used this CA to differentiate different multivalent V_HH constructs and their solution behavior.

Differentiation between multivalent V_HH constructs with varying linker lengths

Candidate selection between two multivalent V_HH molecules (Figure 3a, b; highlighted by a black arrow), that have identical domain topology, but that have linkers that vary in length is not straightforward. Linker length can modulate potential interactions between individual domains but when these interactions are absent, each domain can act independently and behave as if they were isolated in solution.^58,59 Importantly, these dynamic events can influence peak parameters within NMR data. The large variance in the magnitude of the NMR peak parameters for multivalent V_HH domains (Figures 2, 3) is a strong indicator that their dynamics are complex. For multivalent V_HH constructs, interdomain interactions are not desired because they can affect target binding and solution stability.⁵⁶ Therefore, it would be advantageous to identify multivalent V_HH domains that do not cluster near their individual V_HH counterparts and be able to flag them for the existence of interdomain fluctuations. In Figure 3a and b, the two molecules are labeled as short and long. Both molecules maintain identical domain topology, but the long molecule contains a four-fold longer linker as compared to the short molecule. Earlier studies (Figure S3) revealed that the short linker construct sampled states that were larger in the effective MW. This solution behavior was devoid in the longer linker construct (Figure S3), and is attributed to intermolecular interactions that restrict individual domains to tumble independently, causing it to sample configurations that are larger in apparent size. This solution-based characteristic was flagged as a risk for the short linker construct. These molecules were included in order to pressure test the sensitivity of the CA approach. Indeed, the short linker along the PC1/PC2 trajectory lies toward the right of the long linker (Figure 3a), recapitulating the same observation as was seen from orthogonal datasets (Figure S3). An advantage of the CA method is that the peak parameters can be sensitive to lowly populated (<10%) transiently sampled self-associated configurations that interconvert (fast kinetic transitions; microseconds-milliseconds).^42,60 However, establishing increased confidence in the ability to differentiate these possible features with the CA method are needed. To that end, we further investigated the higher order PCs in order to evaluate if additional confidence regarding dynamic features of molecules can be differentiated.

The first four PCs accounted for 90% of all input data. Although, PC1 and PC2 are strongly correlated and report on the shape and size of different Ab(L/F)s (Figure 2), we investigated if the third and fourth PCs provide additional differentiation toward solution behavior. In Figure 3b, PC3 and PC4 are plotted against each other and are significantly less correlated compared to PC1 and PC2. This indicates that these higher-order PCs may report on different aspects of the peak parameters that were input for this CA. Visual inspection of Figure 3b shows a regrouping as to how the different classes of Ab(L/F)s cluster. Fab and multidomain V_HH molecules fall to the right of the plot, single domain V_HH molecules cluster to decreased values and mAbs group to having no significant variance within PC3 and PC4. Since the input data are sensitive to the amplitude and timescale of motional processes, we hypothesized that the higher order PCs may reflect features that are encoded from dynamic events that add to LW_app rather than size and shape information.

We tested whether higher order PCs may reveal dynamic features by simulating different distributions of synthetic methyl relaxation rates and relative peak intensities and then processing the data through our numerical CA (see the Supplemental Figure S4). In brief, distributions of relaxation rates of ¹H and ¹³C nuclei were randomly generated for three different classes of Ab(L/F)s including individual V_HH domains, Fabs, and mAb molecules. A bimodal distribution was used to model the relaxation behavior of molecules with increased internal or interdomain dynamics. The simulation parameters we used were rudimentary and do not completely capture the intricate complexities that can modulate Ab(L/F)s (Figure S4). However, the simulation did recapitulate a critical observation that was made from the experimental data. A positive correlation between PC1 and PC2 (Figure S4; r = 0.86) from the simulated dataset was recovered and is very similar to the experimental reference datasets in Figure 3a. This again supports that LW_app is an appropriate metric to capture size/shape differences between Ab(L/F)s. Projection of PC3 against PC4 from the simulated data revealed plots with increased scatter for the larger Fab and mAb molecules that were modeled to maintain dynamic processes and fell to the right of the PC3/PC4 plot. The short and long linker multivalent (Figure 3a, b) V_HH domains are also further separated along PC3 and PC4, providing additional confidence that this protocol can discern interdomain interactions. The absence of features, within PC3 and PC4, for the molecules tested here, may indicate that dynamic events are not well sampled within their natural abundance [¹H, ¹³C]-HMQC datasets. Interestingly, Fab molecules that were included in the reference set display significant values for PC3 and PC4 (Figure 3b), indicating that some dynamic process may affect these select Fabs. Ordering by decreasing magnitude of PC3 reiterates this observation, providing a way to rapidly identify differences between the different class of Ab(L/F)s (Figure 3c). Rank ordering highlights outliers that differ amongst a single class of Ab(L/F)s. In particular, the long chain multivalent V_HH domain (Figure 3c; blue point with an asterisk) has diminished dynamics as compared to the shorter chain multivalent V_HH constructs, which may suffer from aberrant interdomain interactions. A combination of rank ordering PC3 by magnitude and the higher order components PC3/PC4 provides increased confidence in defining solution behavior for multivalent V_HH domains with differing linker lengths.

A random sampling of NMR peaks for covariance analysis is sensitive to change in formulation

Development of Ab(L/F)s requires assessments of protein behavior in a variety of different buffers. Techniques that detect changes of a molecule, under different formulated conditions are required in order to identify the best means to develop a molecule. 1D NMR-based methods have been previously used to measure changes to a mAb because of formulation.³¹ The PROFILE method was applied to evaluate one mAb, under two unspecified conditions, and the authors reported that small changes in a correlated parameters (4%) could be discerned from 1D ¹H NMR spectra. We propose that 2D NMR fingerprints processed through this approach are also sensitive to certain formulation conditions because differences in the NMR data are related to physical mechanisms that depend directly on shape/size, reorientation time or molecular dynamics (Figure 1, S1; Figure 2). We evaluated three samples of a mAb (called Test-mAb) formulated under three different conditions labeled A (Figure 3d), B (Figure 3e) and C (figure 3f), but with the mAb concentration kept constant. Each 2D dataset was processed through this CA, compared alongside the original reference set of Ab(L/F)s, and then the PC1/PC2 trajectory was plotted. Each new dataset was added as an additional input array. Importantly, each condition did not induce any significant change to the rank of the other Ab(L/F)s across PC1 and PC2 (Figure S5). This indicates that the addition of one new molecule still preserved the PC1/PC2 trajectory, permitting quantitation of any new mAb. It is possible to add all three conditions to the input array simultaneously for the calculation, but we prefer to add one at a time because their behavior is not known a priori and comparison of the single dataset is easier. Each condition for the Test-mAb plotted in a unique position along PC1/PC2 (Figure 3d-f; enlarged black point), with conditions A to C falling in ascending order along the PC1/PC2 trajectory. This indicates that the three different conditions perturb Test-mAb in such a way that their resonance peak properties are affected to different amounts. Importantly, the magnitude of PC1/PC2 eigenvectors reported changes greater than 50% between conditions A through C. The CA approach revealed that condition A causes the Test-mAb to have an apparent behavior similar to molecules that are smaller than other mAbs, while condition C causes the Test-mAb to display characteristics comparable to the mAbs included in the reference set. It should be stressed that the conditions tested here likely stem from changes in the molecule’s overall tumbling time or changes in dynamic processes that can contribute to LW_app. We propose that this level of differentiation is an important feature because it allows a molecule’s behavior to be flagged, due to formulation. For example, if a molecule falls along an extreme of the PC1/PC2, the magnitude of the eigenvectors can be recorded for each formulation and used for comparison. Also, it could be useful to use this approach to flag molecules that have extreme deviations on PC1/PC2 from their known MW_eff. This is because values that are too large (e.g., aggregation, slowed behavior) or too small (e.g., degradation) could be due to aberrant effects. This could also be accentuated by the effect of excipients that alter viscosity of the solution. Importantly, a change of multimeric or aggregated state could be identified by this approach, as that is highly sensitive to the input LW_app.

Discussion

Natural abundance NMR of Ab(L/F)s have a growing role in the qualification of different molecules.^30,41 This approach provides the ability to monitor the behavior of Ab(L/F)s at atomic resolution, under near-native or different formulated conditions. We have demonstrated that a collection of [¹H, ¹³C]-HMQC datasets, of different antibody and antibody-like molecules (Figure 1) can be directly compared. A numerical approach that uses resampling of resonance peak parameters coupled with PCA creates a common basis that projects and simplifies these datasets. An advantage is that this approach is sensitive to the molecular size, shape and dynamics of the input molecules. Random sampling of resonance peak perturbations through this CA method identified an excellent correlation along the PC1/PC2 trajectory (Figure 2) that is well correlated with the relative shape/size of different Ab(L/F)s (Figure 2). Analysis of other PCs also revealed some sensitivity to differences to interdomain interactions between multivalent V_HH molecules (Figure 2a-c). A comparison of the first four modes highlighted that this approach can flag undesired behavior in multivalent V_HH constructs (Figure 3). Simple simulations of LW_app qualitatively corroborated these findings (Figure S4), but further investigations are required to reveal the exact nature that define these higher order PCs. We also evaluated this method’s ability to detect changes of the same molecule formulated under different conditions (Figure 3d-f). The formulations that were tested induced significant changes in LW_app and caused the Test-mAb (Figure 3d-f) to project in different locations along the PC1/PC2 trajectory. This presents an opportunity to inform on how formulation can affect Ab(L/F)s using data that is sensitive to rotational tumbling properties in solution as opposed to more traditional techniques (e.g., laser-based, UV and/or optical).

The selection of molecules that were used here establish a proof of concept on how parameters from 2D NMR datasets can be quantified across different Ab(L/F)s. However, this technique could be enhanced with a larger number of unique but qualified molecules within the reference set. The addition of molecules with different behavior will enhance the precision to which new molecules, of unknown behavior, can be mapped. A limitation of this approach is that only processes with a significant amplitude can modulate LW_app and the relative peak intensity. For example, processes like self-association, multi-domain protein fluctuations, or intrinsic large amplitude changes like slowed loop dynamics could contribute to these NMR observables, but if these processes occur on a fast timescale (<μs) or too slow (>100s ms), they will not have significant impact on the input peak parameters. We anticipate that more accurate modeling of Ab(L/F)s through coarse-grained or atomistic molecular dynamics simulations, coupled with back calculation of NMR-based observables, could reveal the sensitivity of this CA to the intrinsic dynamic properties of Ab(L/F)s.

Our method requires protein concentrations greater than 1 mg/ml for small V_HH domains and greater than approximately 20 mg/ml for mAbs (Table S1), although high-quality natural abundance datasets could be recorded with the NISTmAb at only 10 mg/mL. Based on these sample requirements, this technique could have applications in late discovery/early process development, whereby it could be plugged into existing workflows or used to obtain orthogonal insight into a molecule’s behavior (Figure S6). The sensitivity to subtle differences and weak physical processes makes it attractive when chromatographic or light scattering type experiments cannot discern differences. Another advantage of NMR-based solution characterization of Ab(L/F)s is that the sample requires no modification to its environment, unlike chromatographic techniques, which subject the Ab(L/F) to a column matrix environment that can induce unwanted interactions with the matrix. This presents a unique opportunity for our approach to be applied to coformulations where resonance parameters could be monitored in bulk (all resonance peaks) or for each molecule independently (assuming a sufficient number of unique resonances are observed; see Materials and Methods). That would allow each molecule’s behavior, under a given formulation, to be compared and quantitated. Bulk NMR relaxation rates can be measured via integration of different spectral regions for an Ab(L/F)s spectrum. For molecules that exhibit uniform relaxation behavior, it is possible to draw conclusions because the rate measured is an accurate average value for the molecule. However, Ab(L/F)s are dynamic in nature and the use of a bulk relaxation rate fails to capture an accurate assessment of their molecular behavior.³⁵

Comparing multiple Ab(L/F) types to one another with this parameter alone is not trivial. 1D PROFILE and Diffusion NMR are rapid tools that can be used to obtain a general idea about sample properties (e.g., folded-ness and size),^31,34 but the data derives from 1D information that can cause complications in analysis because of signal overlap. Here, we aimed to expand upon these efforts by taking advantage of high-quality 2D methyl spectra that are accessible for a broad range of Ab(L/F)s (Figure 1a-d). First, use of the NMR peak parameters from these datasets ensures that the input data are collected at residue-level resolution. Then, use of the CA method presented here allows a foundation to be established regarding how a quantitative comparison can be made between different 2D methyl spectra. Typically, using high-resolution NMR data for decision-making when comparing non-identical molecules or when trying to discern differences in the context of formulation has been challenging. We propose that our approach can alleviate these difficulties by providing an orthogonal platform for quantitative comparison of Ab(L/F)s that can enable decisions around the discovery and developability of these types of molecules. Novel scaffolds, including non-antibody molecules (e.g., Alphabody,⁶¹ DARPin,⁶² biased cytokines^63,64), are attractive options for biological signaling modulation. Understanding their behavior will require additional spectroscopic tools that can illuminate their differences as compared to more conventional scaffolds. A random sampling of NMR Peaks for covariance analysis is a tool that functions independent of scaffold. We envision that this tool can contribute to the decision-making around the discovery of molecules that may have therapeutic potential.

Materials and methods

Nuclear Magnetic Resonance Spectroscopy: All antibody and antibody-like molecules were queried with 2D SOFAST-HMQC experiments. A total of 50 complex points were used in the indirect (t_1,max = 13.3 ms) with 1,024 complex points for the direct dimension (t_2,max = 63.9 ms). For each point in the indirect dimension, 2,048 transients were acquired. All 2D spectra were collected at 300 K on a spectrometer operating at a Larmor frequency of 800 MHz. All NMR data were processed using the SMILE algorithm within NMRPipe.⁴⁷ Although, the SMILE algorithm is used for non-uniformly sampled data matrices, this approach still improves the quality of processed spectra on a fully sampled data matrix. Prior to application of the SMILE algorithm, the indirect dimension was forward predicted by a factor of two. A shifted sine-bell was applied in both the direct and indirect dimensions and spectra were zero-filled to 4k points in both dimensions. This was done to improve the detection of peak amplitudes by the auto peak picker. All of these operations were performed in the NMRPipe software suite. Data was acquired for the variety of samples (Figure S1) that maintained concentrations between 1 and 25 mg/mL. All spectra were collected on the same NMR instrument, acquisition parameters and were processed identically. For data collected across different static magnetic field strengths (e.g., 1 GHz, 800 MHz, 600 MHz) a scaling factor could be applied to facilitate their comparison. However, we suggest that data be collected across the same static magnetic field strength. Small differences in base frequency will not translate to significant differences in LW_app and will not impede application of this CA method. Acquisition parameters like the total number of points in the direct dimension and processing parameters should be preserved. Intensities were normalized when each spectrum was initially peak picked, therefore small variations in the number of transients per point collected will not affect downstream analysis. Resonances were picked from within the 2D NMR datasets using an automated picking protocol that used signal-to-noise ratio constraints in order to distinguish resonances from spurious noise or buffer peaks. Peaks were required to be at least two times the baseline noise of the spectrum. A constant spectral window that spans the methyl region typical for proteins samples was used for all spectra analysis. For each spectrum, the peaks were picked within a spectral range of −0.5 to 1.9 ppm (8 ppm to 31 ppm) for the ¹H (¹³C) dimensions. For each peak that was detected within a given spectrum, the LW_app in both ¹H and ¹³C dimensions and the relative peak amplitude (I_peak/I_max; where I_max denotes the peak within the spectral window with the largest peak amplitude) were kept for input into the resampled PCA protocol. Peaks arising from buffer components that overlapped within the constant spectral window were manually inspected and removed from the respective spectrum’s final peak list.

Computational Analysis: The calculation of PCs used a random subset of 60 peaks per calculation. The input arrays contained the ¹H, ¹³C linewidths and relative intensities for the 60 peaks. Note that even for the lowest MW Ab(L/F)s used here (individual V_HH domains), at least 50 methyl resonances alone are expected. This means that some methylene peaks may have been detected within this spectral window. The total number of peaks per spectrum that were chosen at random was determined from the Ab(L/F)s spectrum with the least number of resonances (individual V_HH). The number of peaks selected per iteration was set to approximately 70% the total number of observable resonances. We also evaluated if reducing the number of peaks, selected at random, would affect the results derived from the calculation. We evaluated a range of peak subset selections from 30% to 70%. Even reducing the minimum subset to 30% of the observable resonances was well correlated with the selection set at 70% (PC1; r = 0.99, PC2; r = 0.8, PC3; r = 0.77, PC4; r = 0.86). A key component of the calculation is ensuring that a sufficient number of selections is made in order to minimize the effect of individual outlier LW_app values. For each collection of NMR spectra, PCs were calculated over 10,000 resampled input peak lists. Although, the relaxation behavior of methylene groups is different from methyl groups, we found that a small number of these methylene resonances have small significance during the calculation as the peak selection is repeated many times and is random. Therefore, the effect is averaged out within the distribution of eigenvectors that are solved during the resampling of PCs.

Individual methyl resonances can overlap in larger Fab and mAb molecules. Here, we choose to keep possibly overlapped peaks in the calculation in order maximize the number of possible input peak parameters. When the overlap is very severe (>2–3 methyl resonances), an overlapped peak can give an exaggerated LW_app and can generate an outlier LW_app. We found that resampling the data 10,000 times was sufficient to remove their impact on the distribution of extracted eigenvectors. The first four PCs captured approximately 90% of all the input data (Figure S2) and the magnitude of the eigenvalues were well preserved throughout all iterations of the calculation. The distribution of eigenvectors from each sample and for each PC is fitted to a gaussian function (Figure S6). This value is used for the variance of that sample within that PC. Here, we are concerned with the magnitude of the eigenvectors and not their eigenvalues. During the PCA calculation, the covariance matrix is maximized. The output eigenvectors from each individual PCA calculation can have altered signs depending on the iteration of the calculation. However, the sign only induces a change in the direction, which is arbitrary for our purposes here, as we are interested in the variance, which stays the same (i.e., the magnitude of the eigenvectors is preserved). Therefore, here only the absolute value of the eigenvectors was retained for further analysis. All calculations were performed with in-house written Python scripts that use the Numpy/Scipy libraries. A copy of the source code for “Random Sampling of NMR Peaks for Covariance Analysis” is available in the Supporting Information.

The test set of Ab(L/F)s molecules was probed by DOSY/Diffusion NMR following previously published protocols.³⁴ Diffusion NMR was used to determine the translational diffusion coefficient (D_t) for a given molecule. Diffusion NMR experiments were conducted using conventional pseudo 2D Pulse Field Gradient (PFG) filtered 1Ds where the gradient amplitude was varied from 5% to 95% of the maximum value tolerated by the probe and a total of 32 different gradient amplitudes were collected per sample. The delay between recovery gradients was set to 300 ms where each gradient was applied for 1.5 ms. The maximum gradient amplitude was set to 58.9 G⋅cm⁻¹. Each 1D was adopized using a shifted sine-bell and zero-filled prior to Fourier transform.

The Diffusion NMR spectrum was subsequently Laplacian transformed, using the processing protocol in Bruker TopSpin®. The translational diffusion coefficient was determined by plotting all diffusion values across the 1D spectrum and fitting the values to a Gaussian distribution. The average was taken as the translational diffusion coefficient. An offset between the measured translational diffusion coefficients and the real translational diffusion coefficients can exist. This offset was corrected for by measuring a set of standard molecules, with known diffusion coefficients/R_h, plotting them against one another and fitting to a straight line (slope set to 1). An offset of −0.167 [log₁₀(D_t)] was determined and applied to the measured diffusion coefficients. This corrected value was then used to back calculate R_h using the Stokes-Einstein equation. Viscosity was assumed to be constant for all samples studied here. This assumes that the molecule is spherical and was applied uniformly across all Ab(L/F)s studied here. This is a simplified view as different topologies may require shape-specific models. Future incorporation of computational models that describe hydrodynamic properties, from structure, could improve the correlations observed in Figure 1. All diffusion experiments were conducted at the same temperature that the 2D datasets were collected.

Abbreviation

Table

1D	one-dimensional
2D	two-dimensional
Ab	antibody
Ab(L/F)	Antibody-like and Antibody fragments
CA	Covariance Analysis
DOSY	Diffusion Ordered Spectroscopy
Dt	translational diffusion coefficient
Fab	antigen binding fragment
H(S/M)QC	Heteronuclear (Single/Multiple) Quantum Coherence
IgG	Immunoglobulin G
LW	Linewidth
mAb	monoclonal antibody
ms	millisecond
MW	molecular weight
NMR	Nuclear Magnetic Resonance
PC	Principal Component
PCA	Principal Component Analysis
ppm	parts per million
Rh	hydrodynamic radius

Acknowledgments

The authors would like to acknowledge colleagues across Computational and Structural Sciences, Discovery Biologics, Bioprocess, Protein Sciences and Daniel Wyss, Mikhail Reibarkh, David Roush and Juan Alvarez for critical reading of the manuscript.

Disclosure statement

The authors report no conflict of interest.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website