Introduction

There is great interest in the biotechnology field in studying high concentration therapeutic monoclonal antibody (mAb) solutions. Monoclonal antibodies are typically administered by IV or subcutaneous (SC) injection at high concentrations (Rituxan: 10 mg/ml; Avastin: 25 mg/ml; Simponi: 100 mg/ml; Humira: 100 mg/ml; Herceptin: 150 mg/ml; Xolair: 150 mg/ml). This requires that mAbs be highly soluble and stable under formulation conditions. Analytical Ultracentrifugation (AUC) sedimentation velocity (SV) is highly appropriate for studying mAbs because it can reveal the presence of unfolding, dissociation, aggregation and undesirable reversible association over a wide range of concentrations (Berkowitz 2006; Philo 2009; Liu et al. 2005). To address the needs of therapeutic protein development, experiments at high concentration and in complex solution environments are required (Shire et al. 2004; Liu et al. 2015). The fluorescence detection system (FDS) developed by Tom Laue and colleagues (MacGregor et al. 2004; Kroe and Laue 2009; Kingsbury and Laue 2011) and, until recently, available from Aviv, is an extremely useful technique for studying high concentrations and heterogeneous systems like serum. FDS experiments with mAbs in serum can reveal unusual interactions that do not occur in typical formulation buffer conditions (Demeule et al. 2009). Nonideality is represented by the following phenomenological equations for sedimentation and diffusion (Fujita 1962).

$$s= \frac{{s}^{\mathrm{o}}}{1+{k}_{s}c} \,\text{and}$$
(1)
$$D= \frac{{D}^{\mathrm{o}} (1 + 2B{M}_{1}c )}{(1 + {k}_{s}c)}$$
(2)

In these equation so and Do are values extrapolated to zero protein concentration, c is concentration in mg/ml, ks is hydrodynamic nonideality and BM1 is thermodynamic nonideality (the 2nd virial coefficient B times molecular weight), both expressed as ml/mg to match the units of concentration. The antibody field has used orthogonal techniques like static light scattering (SLS) to measure a second virial coefficient, BM1, and DLS to measure concentration dependence of diffusion, kD. This is typically summarized as kD = BM1 − ks (Harding and Johnson 1985) suggesting concentration dependence of diffusion involves both a hydrodynamic and thermodynamic component. Many studies combined light scattering, AUC and DLS to investigate the behavior of antibody solutions, emphasizing the complementary information gathered on nonideality and, where present, association (Solovyova et al. 2001; Saluja et al. 2010; Yadav et al. 2011a, b). The goal in part is to survey formulation conditions for the presence of aggregation (Saluja et al. 2010), the avoidance of high viscosity (Yadav et al. 2012), or conditions that favor crystallization (Solovyova et al. 2001). The discussion often focuses on the algebra, do ks + kD equal BM1, and the problems with subtracting small values, BM1 − ks, that nearly cancel. This algebraic expression shown above is valid only in the absence of association. Thus, if ks, kD and BM1 are negative, a nonideal association model should be implemented, such that the nonideality and the association contributions can be accounted for separately. This is generally not feasible by DLS and/or not implemented by SLS methods. The analysis of weak, nonideal association and the extraction of both nonideality and association constants by AUC methods is challenging; it cannot be done with distribution analysis (Rowe 2011; Wright et al. 2018a). On the other hand, rigorous analysis requires nonlinear least squares (NLLS) global fitting of AUC SV data sets to an explicit model that includes both nonideality (ks and BM1) and association. This approach is feasible in SEDANAL because of the ModelEditor feature that allows one to construct any number of simple or complex nonideal associating fitting models (Stafford and Sherwood 2004; Correia and Stafford 2015; Sherwood and Stafford 2016).

Here we present a series of simulations of the impact of ks and BM1 on high concentration antibody solutions of the therapeutic IgG Simponi. We present this over a range of values (0–100 ml/g) individually, ks and BM1 alone, and together, both ks and BM1, to allow a basis set of features that reveal their impact on SV boundary shape. We then choose realistic ks and BM1 values, and simulate a wide range of protein concentrations (1–120 mg/ml) that overlap Simponi’s SC therapeutic dose, 100 mg/ml, and demonstrate the range of data feasible for experimental determinations. Finally, in three recent studies, we showed that all mAb IgG’s weakly self-associate and hetero-associate with other classes of IgG (Wright et al. 2018a, b; Yang et al. 2018). Thus, we present here simulations of nonideal, weak association to demonstrate the challenge of dissecting out basis sets for ks, BM1 and association. This analysis has been possible with SEDANAL for a decade. The full power of the method required the implementation of features first described by Todd and Haschemeyer (1981), and this only recently was included in SEDANAL version 6.97 and above. These simulations provide graphical and conceptual insight for AUC users, and encourage the use of rigorous and robust global direct boundary fitting methods for nonideal, weakly associating systems.

Methods

Simulations were performed with SEDANAL, version 6.97 or later, a software package developed to fit AUC data for complex, nonideal associating systems (Stafford and Sherwood 2004). The models simulated include ks (hydrodynamic nonideality), BM1 (thermodynamic nonideality), and weak self-association, represented here as dimerization, K2. The mAb modeled is golimumab or Simponi, with a molecular weight of 146,909, an extrapolated so = 6.6 s, an experimentally measured buoyancy (1 − νρ) of 0.26623. Simulations in absorbance mode were done with no added noise. Since in silico absorbance or interference experiments are not limited by optical constraints, we adjusted the extinction coefficient for different path cells and present normalized plots. Data were simulated with 1500 points, and output with 650 points between 5.9 and 7.2 cm to match typical absorbance experiments.

Simulations of FDS tracer mode were done by setting the signal strength (extinction coefficient × quantum yield, referred to here as ε1) for monomer at 1 mg/ml to 1000. To produce a constant signal in tracer mode each apparent extinction coefficient was calculated as 1000/co, where co is in mg/ml. For example, cell 2 is 5 mg/ml so ε2 = 1000/5 = 200; cell 3 is 10 mg/ml so ε3 = 100; cell 9 is 120 mg/ml so ε9 = 8.333. All FDS simulations have added Gaussian noise of 10 corresponding to 1% of signal. SEDANAL global fitting of these data sets uses ~ 60–100 scans per data set, which amounts to up to 200,000 total data points, using the Levenberg–Marquardt algorithm to minimize standard deviation. Experimental in vitro or in serum optical constraints are addressed in the discussion.

Values of ks and BM1 were varied over a wide range (0–100 ml/g) to investigate and display their impact on boundary shapes. Experimentally we have shown that Simponi exhibits ks and BM1 values of ~ 10 ml/g (Wright et al. 2018b). (SEDANAL internally uses concentration units of mg/ml and thus ks and 2BM1 have units of ml/mg to make ksc or 2BM1c dimensionless (Eqs. 1 and 2). However, the field uses ml/g (Rowe 1977; Harding and Johnson 1985; Solovyova et al. 2001; Saluja et al. 2010; Yadav et al. 2011a, b) and thus we plot data vs mg/ml, but report data in units of ml/g (Table 1). This avoids reporting numbers with lots of digits (0.010 ml/mg vs 10.0 ml/g).) Thus, we also separately constrained ks and BM1 to 10 ml/g and varied total concentration from 1 to 120 mg/ml. Weak self- and hetero-association of therapeutic mAbs has been demonstrated (Wright et al. 2018a,b; Yang et al. 2018), and is assumed to be dimerization because of the small extent of reaction. The goal of these simulations was to present graphical basis sets that demonstrate the impact of ks and BM1 on boundary shapes. Thus, simulations are initially presented over a wide range of values to give a picture of what nonideality looks like in an experimental setting. Simulations were also performed at fixed values of ks and BM1 and as a function of concentration (1–120 mg/ml) to mimic therapeutic samples and to establish what is required to extract experimental values of the parameters. We observed that all mAbs self-associate and hetero-associate with other antibodies. To explore this, nonideal association (ks, BM1, K2) simulations as a function of small K2 values are presented. In addition, all mAbs contain small amounts of dimeric and trimeric aggregates, and these are included in the model to mimic realistic heterogeneity.

Table 1 SEDANAL analysis of high concentration FDS data
Full size table

All ks and BM1 values are provided as matrices, with elements, kij and BijMi, representing self- and cross-term interactions (Table S1 and S2; Correia and Stafford 2015; Wright et al. 2018a). This represents a family of phenomenological equations for each component that explicitly includes all species concentrations. For example, s1 (or sA) can be written

$${s}_{1}= \frac{{s}_{1}^{\mathrm{o}}}{1+{k}_{11}{c}_{1}+{k}_{12}{c}_{2}+{k}_{13}{c}_{3}}$$
(3)

where c2 corresponds to dimer concentrations, c3 to trimer concentrations, and k12 and k13 are cross-term nonideality terms reflecting the effect of dimer and trimer concentrations on monomer s1. Other species like reversible dimers or irreversible dimers and trimers have similar expressions (see Supplemental equations S1–S4). In addition, D1 (or DA) can be written

$${\mathrm{D}}_{1} = \frac{{\mathrm{D}}_{1}^{\mathrm{o}} (1+{2B}_{11}{M}_{1}{c}_{1}+{2B}_{12}{M}_{2}{c}_{2}+{2B}_{13}{M}_{3}{c}_{3})}{(1+{\mathrm{k}}_{11}{\mathrm{c}}_{1}+{\mathrm{k}}_{12}{\mathrm{c}}_{2}+{\mathrm{k}}_{13}{\mathrm{c}}_{3} )}$$
(4)

with similar expressions for dimers and trimers (Correia and Stafford 2015; Stafford 2016). Note we assume the cross terms kij and BijMj for dimers and trimers are the same as monomers on a weight scale; this is reasonable since experimentally mAb aggregates have the same frictional ratio, f/fo (Philo 2003), but this is adjustable in the matrix if required. (See Supplemental Methods section for more details.)

SEDANAL simulates and fits SV data by the method of Claverie (1976) which uses finite-element solutions to the Lamm equation, the partial differential equation that describes simultaneous sedimentation and diffusion in a sector shaped cell (Todd and Haschemeyer 1981).

$$(\partial {c}_{i}/\partial t{)}_{r}=-\frac{ \partial }{r\partial r} \left[{c}_{i} {\omega }^{2}{s}_{i}{({c}_{i})r}^{2}-{D}_{i}({c}_{i})r{\left(\frac{\partial {c}_{i}}{\partial r}\right)}_{t}\right]$$
(5)

The si and Di terms are concentration-dependent and defined as above (Eqs. 3 and 4) to reflect nonideal behavior. If mAb samples are heterogeneous, paucidisperse systems and contain small amounts of dimeric or trimeric aggregates, simultaneous Lamm equations for the dimer and trimer are solved (Correia and Stafford 2015). Reversible association is included by adding a relationship that links dimer concentration to monomers through equilibrium or kinetic rate constants (Stafford and Sherwood 2004). Figure 1 presents an example of a SEDANAL fitter model for simulation. All antibody samples contain small amounts of irreversible dimer and trimer aggregates, set here to 10% and 5% weight fraction. In addition, antibodies can undergo reversible nonideal association represented here as monomer–dimer (AA2) plus ks and BM1. Thus, in this model system, AA2 refers to reversible dimer formation and BC or dimer aggregates and trimer aggregates refers to irreversible dimers and trimers.

Fig. 1
figure 1

Screen shot of a SEDANAL AUC SV simulation of Simponi at 40 K rpm, 10 mg/ml, with ks = 0.01 ml/mg, BM1 = 0.01 ml/mg, and K2 = 5000 M−1, with 10% dimer and 5% trimer. Other details are visible in the figure

Full size image

High concentration solutions to the Lamm equation are intrinsically unstable, primarily due to numerical instability in the base region of the cell (SEDANAL detects this as a “Check Grid” error). Todd and Haschemeyer (1981) were the first to develop curve-fitting methods using Claverie’s rapid numerical solutions to the Lamm equation; the Claverie method was for comparing simulations to real data (Claverie 1976). We constrained parameters to deal with this instability. The most significant of these are defined as ε (epsilon) and η (eta), where s/so = 1 − ε and D/Do = 1 − η (Eqs. 3 and 4). For realistic parameters, meaning ks and BM1 equal and positive, D/Do < 2 (Figure S1A); other choices give less restrictive ranges (Figure S1B). These constraints are provided to the SEDANAL fitter to set the upper limit on D/Do. SEDANAL also allows concentrations at the base to be limited.

We restrict ks and BM1 to positive values because negative values imply association, and therefore, nonideal associating systems are modeled by an explicit nonideal, association scheme (see Fig. 1). Otherwise, the magnitude of ks and BM1 reflect the sum of multiple attractive and repulsive interactions and are difficult to interpret (Laue 2012; Wright et al. 2018b; Laue and Shire 2020). Sedimentation and diffusion coefficients vary across the cell due to nonideality and to local changes in concentration as the boundary components sediment and fractionate (Figure S2). This fractionation or boundary separation has been called demixing by Kingsbury and Laue (2011). This variation in the effect of nonideality is especially true in the base region, where concentrations can increase by two or more orders of magnitude, and calculations of concentration become numerically unstable. SEDANAL has been upgraded to handle these effects (Todd and Haschemeyer 1981), and all the data simulated and fit in this study exhibited stable solutions.

Simulated data are analyzed by Wide Distribution Analysis (WDA) an option in SEDANAL, and plotted as s*g(s*) vs log(s*). This transformation follows from the definition of g(s*) = dc/ds; thus s*dc/ds* = dc/d(ln(s*)) (Stafford and Braswell 2004; Sherwood and Stafford 2016). WDA was developed for multi-speed experiments but we have found it useful for rapid analysis of single-speed runs as well. We plot data as normalized s*g(s*) vs log(s*) for convenient reading of the scale and for direct comparison of trends in the distributions. The WDA method allows all data scans to be included, with early scans providing information about large components (out to 1000 S in these simulations), and late scans providing information about small species. One can choose different radial positions to view sedimentation at different extents of resolution (Figure S3). A standard DCDT analysis (Stafford 1992) cannot mimic this degree of resolution because it is limited to a narrow span of data scans, in this case a region where the monomer boundary is near 6.4–6.5 cm. Nonetheless, DCDT+ analysis (Stafford 1992; Philo 2006) can be compared with s*g(s*) (Figure S4) to reveal the precision and advantages of both methods. DCDT+ and WDA generate sw and integrated signal area under the curve from results extrapolated to the initial concentration, co (Kegeles and Gutter 1951; Patel et al. 2018). Both methods are used to construct 1/sw vs c plots (Figs. 6 and 10) to explore the limitations of linear graphical analysis.

Results

Figure 2 shows a family of simulations performed at 10 mg/ml and over the ks range 0–100 ml/g. This corresponds to a k11c1 range, ml/mg times mg/ml, of 0.01 to 1.0 at co, or an initial sapp range of 6.6–3.3 s. Thus, the monomer boundary slows, and, since (1 + k11c1) is also in the denominator for D (Eq. 2), the monomer boundary is sharper. Note, the concentrations vary with radial position and time in the simulation due to demixing and radial dilution, and thus the exact sapp correction is a function of radial position and time throughout the cell and involves contributions from monomer, dimer and trimer species, si = (1 + ki1c1 + ki2c2 + ki3c3) (Figure S2). The boundaries for dimer and trimer sediment in the plateau region of the monomer, and are thus dominated by the kij cross terms for the larger monomer concentration. Note that the aggregate peak widths do not change significantly because they are in the presence of a nearly constant concentration of monomer. This is in contrast to the monomer boundary where low concentration on the centripetal side of the boundary speeds up sedimentation, and high concentration on the centrifugal side slows sedimentation leading to sharpening, increased concentration and a negative concentration gradient. This is the classic Johnston-Ogston (JO) effect (Johnston and Ogston 1946; Correia et al. 1976, 2016). Oddly enough, even at 10% dimer, the JO effect is evident in the monomer boundary as it interacts with the dimer boundary (Figure S2). The JO effect is more pronounced at both high ks and high concentrations since it is due to a product, ksc.

Fig. 2
figure 2

Simulation of ABC or monomer–dimer-trimer model, at 10 mg/ml Simponi plus 10% dimer and 5% trimer irreversible aggregates, where all kij values (self and cross term; Table S1) vary from 0 to 100 ml/g. Data are plotted as normalized s*g(s*) vs s* on a log scale

Full size image

Typical ks values for mAbs are 4–11 ml/g (Wright et al. 2018b; Yang et al. 2018). However, experimentally ks is known to increase dramatically for asymmetric proteins (Creeth and Knight 1965). For example, fibrinogen is 15.5 ml/g; TMV is 33.8 ml/g; myosin is 49–73 ml/g; collagen is 265 ml/g. Pegylated proteins (Li et al. 2012) and IDPs like ELP (Correia et al. 2016) also have very large ks values. Thus, Fig. 2 demonstrates a wide range of expected behaviors for compact to moderately asymmetric proteins.

Figure 3 shows a family of simulations performed at 10 mg/ml and over the BM1 range 0–100 ml/g. Since BM1 only influences D, the boundary spreads, the peak position drifts, but on average does not sediment differently (Sw = 7.12 s ± 0.019 s or 0.27%). This corresponds to a 2B11M1c1 range, ml/mg times mg/ml, of 0.02–2.0 at co, or an initial D/Do range of 1–3 corresponding to the (1 + BijMjcj) term. Note the broadening of the monomer region and an equal degree of broadening of the dimer and trimer regions due to cross term nonideality. This will reduce the resolution of these minor peaks for single sample analysis. Also note this broadening in the monomer region makes the boundary appear more heterogeneous (see “Discussion”).

Fig. 3
figure 3

Simulation of ABC model, at 10 mg/ml Simponi plus 10% dimer and 5% trimer aggregates, where all BijMj values (self and cross term) vary from 0 to 100 ml/g. Data are plotted as normalized s*g(s*) vs s* on a log scale

Full size image

Both ks and BM1 are caused by a combination of excluded volume, charge and shape effects (Tanford 1961; Rowe 1977). As discussed above and in the literature, ks and BM1 have similar magnitudes, and thus this range of values is reasonable in principle (Yadev et al. 2012; Wright et al. 2018b; Chaturvedi et al. 2019; Chaturvedi and Schuck 2019). To demonstrate this similarity we simulate 10 mg/ml Simponi and vary both kij and BijMj equally from 0 to 100 ml/g (Fig. 4). This range of parameter values should impact s values in a manner similar to Fig. 2, but D values now have changes to both the numerator and denominator that offset each other (Eq. 4). In this case the ratio (1 + B11Mc1)/(1 + k11c1) will vary from 1 to 1.5, but as described above the concentrations of all species vary with radial position and time in the run due to demixing and radial dilution, and thus the exact D/Do value is a function of radius and time throughout the cell (Eq. 4). This increase in D/Do causes the dimer and trimer regions to broaden and thus reduces the resolution of the dimer and trimer peaks. Note also that these data exhibit a JO effect at the monomer–dimer interface in a concentration-dependent manner according to the value of kij cross terms (Figure S2).

Fig. 4
figure 4

Simulation of ABC model, at 10 mg/ml Simponi plus 10% dimer and 5% trimer aggregates, where all kij and BijMj values (self and cross term) vary from 0 to 100 ml/g. Data are plotted as normalized s*g(s*) vs s* on a log scale

Full size image

These first three sets of simulations (Figs. 2, 3 and 4) vary nonideality at a fixed protein concentration. Figure 5 shows concentration-dependent simulations of Simponi from 1 to 120 mg/ml plus 10% dimer and 5% trimer aggregates. (Figure S6A presents the DCDT+ g(s*) vs s* version of this plot). Nonideality parameters, kij and BijMj, (self and cross term) are constrained to experimentally reasonable values of 10 ml/g (Wright et al. 2018a). The increase in hydrodynamic nonideality reduces sapp in a manner similar to Fig. 2. The ratio of thermodynamic to hydrodynamic terms (1 + 2BijMjc)/(1 + kijc) in the plateau region of the simulations increases D/Do from 1.01 to 1.58. This is most evident in the dimer and trimer zones where discrete peaks are broadened and more difficult to resolve. As the protein sediments the concentration dramatically increases at the base (Figure S5A), causing significant spreading or back diffusion in the base region of nonideal samples that precludes using radial values too close to the base for WDA (Figure S5B). (It is worth noting that filling a cell to only 6.1 or 6.2 cm reduces the usable radial range for analysis and may be problematic for these highly nonideal systems, especially in FDS experiments where the base is already cutoff optically.)

Fig. 5
figure 5

Simulation of ABC model, at 1 to 120 mg/ml Simponi concentrations plus 10% dimer and 5% trimer aggregates. Nonideality parameters kij and BijMj values (self and cross term) are constrained to 10 ml/g. Data are plotted as normalized s*g(s*) vs s* on a log scale

Full size image

These simulations emphasize three features of highly, nonideal systems due to either large kij values or high concentrations. First, the boundary will be hyper-sharp during sedimentation, due to concentration dependence of sapp, and second, the base region will be very broad due to back diffusion into the plateau region caused by an increase in Dapp. This second feature is entirely due to thermodynamic nonideality and reduces the radial region of data that can be analyzed by typical software packages. Thirdly, minor aggregates typically found in mAb solutions will experience significant broadening making their determination challenging for single sample analysis. In more heterogeneous solutions like serum, the high concentrations at the base cause density gradients that lead to banding of some lipoprotein species. (Analysis of serum samples by AUC FDS methods will be discussed in a future paper). It is worth noting that Figs. 4 and 5 look identical in normalized mode, but Fig. 5 is a concentration series, and thus the signal actually increases dramatically unless done in tracer mode (see below).

AUC data collected as a function of concentration are typically analyzed by linear plots of s vs c or 1/s vs c to extract ks. The first plot is based upon the equation s = so(1 − k′sc), a Taylor expansion of Eq. (1), while the second plot is the linear rearrangement of Eq. (1); 1/s = 1/so + (ks/so)*c. We have previously discussed the relationship between these two coefficients, ks and k′s (Wright et al. 2018b). A 1/sw vs c plot of the data in Fig. 5 is presented in Fig. 6 using sw of the full boundary, the WDA monomer peak positions or DCDT+ sw of the monomer region. The data from WDA and DCDT+ are plotted vs total concentrations or concentrations corrected for radial dilutions (Patel et al. 2018). Matching the sw with the appropriate concentration is the challenge in the analysis of mixtures by a graphical method. The extrapolated so for monomers are reasonably close to the simulated value 6.6 S. The ks values approach the expected value of 10 ml/g but cannot precisely capture the more dynamic nonideality effect of radial dilution and demixing of the species (see “Discussion”). A simulation of just monomers (0% dimers, trimers) gives better but also not perfect estimates of ks (data not shown). The sw data for the full boundary also approach the correct s0w and ks values, but fail to capture the dynamics of nonideality and demixing (Kegeles and Gutter 1951; Patel et al. 2018).

Fig. 6
figure 6

1/s vs c plots of the data in Fig. 5 analyzed with DCDT+ or WDA. The heavy dashed line represents monomer with an intercept at 6.6 s and ks = 10 ml/g. The heavy dot-dash line corresponds to a mixture of monomer–dimer-trimer with a sow = 7.174 and ks = 10 ml/g

Full size image

Choosing the appropriate concentration and sw seems to be the challenge in this linear analysis (Patel et al. 2018). To rigorously deal with sedimentation, radial dilution and species demixing (Kingsbury and Laue 2011), the data should be globally fit to an appropriate model as a function of concentration and time using finite-element solutions to the Lamm equation. To demonstrate the need for global fitting, we simulated FDS tracer data for this same ks, BM1 model over the concentration range of 1–120 mg/ml. A global SEDANAL fit of the data is presented in Fig. 7. The best fit returns s = 6.5994, ks = 10.0 ml/g and BM1 = 9.99 ml/g with an rms = 14.153 which is the expected value for the added noise, 10, increased by the √2 for the ΔC method used by SEDANAL. All the dimer and trimer weight ratios are correct within 1%. Our conclusion is that a full concentration series, analyzed by global direct boundary fitting with a proper ks, BM1 heterogeneous model is the most rigorous and accurate way to analyze high concentration AUC SV data. A bootstrap analysis of these data reveals very tight error bars as expected for simulated data (with 217,450 points; Table 1), but the error bars for BM1 are ~ 6 × larger than for ks. This is consistent with our experimental observation that BM1 is much more difficult to measure than ks by SV analysis.

Fig. 7
figure 7

A global SEDANAL fit of nine FDS data sets to an ABC-ks-BM1 model with 10% dimer and 5% trimer aggregates. Data were simulated with a tracer amount of labeled Simponi corresponding to 1000 counts of monomer, and increasing amounts of unlabeled mAb. To produce a constant signal in tracer mode each apparent extinction coefficient is set to 1000/co where co is in mg/ml. In this in silico FDS tracer mode, extinction coefficient is held constant and concentration is fit

Full size image

We observed that all mAbs exhibit weak self-association properties (Wright et al. 2018a, b; Yang et al. 2018). To include this in the analysis we have simulated a nonideal self-associating model at 10 mg/ml and increasing values of K2 (Fig. 8). At 10 mg/ml and no association (K2 = 0) the curve clearly exhibits nonideality relative to the known monomer s value. At 1000 M−1 the main peak is still slightly nonideal meaning the monomer peak runs at less than 6.6 s. The impact of this is that ks,app will be smaller than expected because it is masked by the association (Wright et al. 2018a; Yang et al. 2018). As K2 increases the main peak shifts to the right indicating association, although it is worth noting the boundary is still nonideal, but now the nonideality is partially masked by weak association.

Fig. 8
figure 8

Simulation of nonideal self-association models as a function of K2 at 10 mg/ml, ks = BM1 = 10 ml/g. Results for K2 values of 0, 1000, 5000, 1e4, and 5e4 M−1 are presented. (Note these K2 values in M−1 correspond to Kd of 0, 73.4, 14.7, 7.3 and 1.47 mg/ml, or K2 (mg/ml)−1 = 2 × K2 (M)−1/146,909) The monomer s value is indicated by the vertical dashed line

Full size image

To investigate this further we simulated a nonideal concentration series for all four K2 values (Fig. 9 presents data for K2 = 5000 M−1) and then plotted all the data as 1/sw vs c, corrected for radial dilution (Fig. 10). Up to ~ 10 mg/ml the data shifts to increasing s values, indicating association, and then shifts to lower s values reflecting an increasing effect of nonideality (Stafford 1980). For weak association, the linear portions of the curve give a reduced apparent ks value consistent with masking of nonideality by association (Wright et al. 2018a). Above 20 mg/ml nonideality becomes dominant, and with increasing K2 values, the slope and intercept approach but do not equal ks and so values expected for the dimer. Nonideality and association are concentration-dependent and thus give different apparent answers at different concentrations. These results support the overlapping impact of radial dilution and fractionation during sedimentation on nonideality and association and the challenge of extracting these values by linear graphical methods. As stressed above, the rigorous analysis of these data requires NLLS global fitting to the primary SV data.

Fig. 9
figure 9

Simulation of nonideal self-association (ks = BM1 = 10 ml/g, K2 = 5000 M−1) as a function of concentration (1–120 mg/ml). The data are plotted as Norm s*g(s*) vs s* on a log scale

Full size image
Fig. 10
figure 10

A plot of 1/sw vs concentration for ks, BM1, K2 data in Fig. 9 and repeated for data simulated with K2 = 1000, 1e4 and 5e4 M−1. The data were analyzed by integration of the full boundaries with DCDT+. The linear portions of the curves were fit to a straight line to extract apparent ks values and extrapolated so values. The dotted line is pure monomer with so = 6.6 S and ks = 10 ml/g. The dashed line is pure dimer with so = 9.9 S and ks = 10 ml/g

Full size image

To demonstrate the requirement for NLLS fitting, we simulated FDS tracer data for this same ks, BM1, K2 model over the concentration range of 1–120 mg/ml for a nonideal system with K2 = 5000 M−1 (Fig. 11). The best global fit of these data returns s = 6.5989, ks = 10.00 ml/g and BM1 = 10.06 ml/g, K2 = 5008 M−1 and an rms = 14.106. The weight ratio of dimer and trimer aggregates are returned correctly to less than 1%. Bootstrap analysis reveals a 3 × larger uncertainty on BM1 vs ks consistent with the observation above (Table 1; 184,340 points). The BM1 and K2 data have a surprisingly small correlation coefficient (R = 0.34). Thus, we conclude that rigorous analysis of high concentration, nonideal, weak, associating systems requires global direct boundary fitting. At present only SEDANAL does global direct boundary analysis of AUC SV data to a nonideal associating model over this concentration range. There are other AUC approaches to analysis of high concentration SV data (Chaturvedi et al. 2018, 2019; Chaturvedi and Schuck 2019) that use a nonideal version of the c(s) method cNI(s). These methods estimate ks by a plot of sw versus concentration, a problematic approach as discussed above, and do not do global fitting. Combining sedimentation equilibrium analysis to independently measure BM1 is certainly a reasonable approach, but it is not as useful for heterogeneous solutions like serum (see “Discussion”). Examples of additional fitting requirements for real experimental data will be presented in a separate publication.

Fig. 11
figure 11

A global SEDANAL fit of nine FDS data sets to an AA2BC-ks-BM1-K2 model with 10% dimer and 5% trimer aggregates and K2 = 5000 M−1. Data were simulated with a tracer amount of labeled Simponi corresponding to 1000 counts of monomer, and increasing amounts of unlabeled mAb. To produce a constant signal in tracer mode each apparent extinction coefficient is set to 1000/co where co is in mg/ml

Full size image

Discussion

The goal of this work is to present graphical distributions from SV simulations of nonideal mAbs solutions at high concentrations up to and exceeding the therapeutic concentration. This allows extraction of features that are direct consequence of hydrodynamic and thermodynamic nonideality. This should provide the opportunity to see the impact of ks and BM1 in SV data sets. Is the boundary hypersharp and running more slowly than expected? Is the base region much broader than typically observed? Is there evidence of negative gradients in the plateau region because of the JO effect, consistent with the presence of mixtures? The phenomenological or empirical parameters ks and BM1 are best determined from their concentration-dependence. It is imperative that experiments be done as a function of loading concentration, the wider the range the better. We simulated data up to 120 mg/ml to exceed the therapeutic dose of Simponi. We also explored the range of 1–100 ml/g for ks and BM1. As discussed above, this range corresponds to compact, globular and moderately asymmetric proteins.

mAbs appear to have ks values in the range of 4–11 ml/g (Wright et al. 2018a, b; Yang et al. 2018), but this is sensitive to conditions, especially salt concentration, and methods of analysis. Weak self-association will mask ks and BM1 values, and thus proper models that include both nonideality and association are required (Wright et al. 2018a). Strongly associating mAbs are rare, but there are reports of much larger ks and BM1 values that may reflect asymmetric complex formation, linear chains of mAbs in an extended conformation (Liu et al. 1995; Hopkins et al. 2018). This discussion is based upon the assumptions outlined by Rowe (1977) that ks is due to swollen volume Vs and f/fo effects according to the equation ks = 2v(Vs/v + (f/fo)3) where v is vbar. The effects of pH and salt on ks are also of great interest. Yang et al. (2018) showed that increasing salt suppresses mAb association and allows ks values to increase and approach values of 10 ml/g. Recently Chaturvedi et al. (2019) and Connolly et al. (2012) reported ks values > 20 ml/g in low salt conditions. Many investigations vary salt concentration and observe a suppression of association and reduction in viscosity (Yadev et al. 2012). Changes in pH have a surprisingly small impact. Direct charge measurements by membrane-confined electrophoresis (Moody and Shepard 2004) show a rather constant mAb charge as a function of pH. They speculate that mAbs buffer charge due to Cl− or anion binding (Yadav et al. 2012; Yang et al. 2018).

The simulations presented here were done in absorbance mode in the absence of added noise. Interference could have worked equally well because in silico has no optical limitations. In practice absorbance experiments on typical mAbs at 280 nm, assuming a generous useable range of 2 OD, are limited to ~ 1.1 mg/ml in a 12 mm cell, to ~ 13 mg/ml in a 1 mm cell. Collecting data at the UV minimum near 254 nm will expand this range, depending upon Trp/Tyr content, by a factor of 2. Interference has a much wider dynamic range and can be useful for AUC SV experiments up to 50 mg/ml (Chaturvedi et al. 2018). Nanolytics has produced a new interference system (Schilling and Krause 2015) that claims to work up to 150 mg/ml. Both of these claims about interference optics are dependent upon accurate optical focusing and conditions, especially speed, since fringe resolution is strongly dependent upon the gradient steepness. As shown in Figure S3, resolution of sedimenting species is strongly speed-dependent. Conditions also imply proper dialysis, the use of meniscus matching centerpieces, and may be influenced by formulation buffers. It is worth pointing out that both absorbance and interference are mostly limited to single component systems like pure mAbs.

Our approach to high concentration experiments uses FDS in a tracer mode, referred to as BOLTS or Biological On-Line Tracer Sedimentation by Laue (2009), where an Alexa-488 labeled mAb is run in a high concentration background of unlabeled mAb (Wright et al. 2018a). This allows no significant limit to the concentration ranges explored, besides solubility, and has the further advantage of being useful in serum, cell, and tissue extracts. As described above, this approach also allows a constant fluorescence signal to be used and thus provides better signal/noise characteristics for data analysis (Husain et al. 2015; Lyons et al. 2013). Our preference is ~ 500–1000 counts (out of 4096 max), which is a function of label efficiency, but generally means 100–200 nM in typical mAb samples. This is the range of FDS data we simulate and fit in Figs. 7 and 11. The results are remarkably good with excellent accuracy and precision. Not all mAb are available to typical users over a range of 1–120 mg/ml. Thus experiments may only be possible over 1–40 or 1–20 mg/ml ranges. Fitting of these more limited concentration ranges also gives excellent results with the expected slightly wider uncertainty limits, especially in BM1 (Table 1). These considerations will be explored in subsequent experimental work.

It has been standard practice in the AUC field to extract ks values from concentration dependent SV analysis mostly through linear plots (Figs. 6 and 10; Kegeles and Gutter 1951; Creeth and Knight 1965; Rowe 1977; Li et al. 2012; Patel et al. 2018; Wright et al. 2018a, b). Although it has been proposed (Solovyova et al. 2001; Wright et al. 2018a, b), it is less common to extract BM1 values from SV analysis. Thus, it is worth asking, what is the consequence of not including BM1 in a direct boundary fit of data like Figs. 7 and 11? To investigate this we refit the AA2BC data in Table 1 without BM1. The 1–120 mg/ml data set returned an increase in rms (rms/rmso where rmso is the fit with BM1) of 23%, while for the 1–40 mg/ml data the rms/rmso increased by 3.5%, and for the 1–20 mg/ml data the rms increased by only 0.8%. To achieve these fits the so values increased slightly while K2 decreased (not shown). The surprising result was that the aggregated dimer to monomer fraction, B/A in the model (Fig. 1), increased dramatically while the monomer fraction decreased. As mention above for Fig. 3, the impact of BM1 is to broaden the boundary making it appear heterogeneous. Thus, in the absence of BM1 the best fit now increases B/A, the aggregated dimer/monomer fraction, to match the apparent heterogeneity of the boundary shape. Thus, not including BM1 in the fit appears to make aggregation seem to increase. To investigate the coupling between BM1 and the other parameters, we repeated the simulations without aggregated dimers and trimers, i.e. without species B or species C, and then compared fitting with and without BM1 (Table 2). Now the rms/rmso values were increased by 5.5%, 19% and 90%, respectively, for the three concentration ranges. This was also matched by slightly larger so values and smaller K2 values, with so and K2 values highly correlated, R = 0.85. Thus, in the absence of aggregated dimers, the omission of BM1 in the fit has a larger impact on the best NLLS fit. These large rms/rmso deviations are consistent with weak correlation between BM1 and the other parameters in the full fits (K2 vs BM1, R = 0.12; s vs BM1, R = 0.18; ks vs BM1, R = 0.14). Since BM1 seems to be nearly orthogonal to the other parameters, s, ks or K2 cannot compensate for the absence of BM1 in the fit. In practice, this will be complicated by the actual presence of molecular heterogeneity due, for example, to variable glycosylation states or mixtures of other mAb conformations.

Table 2 SEDANAL analysis of AA2 ks BM1 K2 = 5000 FDS model
Full size table

As previously reported, the FDS dynamic range in our hands is low nM to low µM. Linearity of signal is the main concern (Lyons et al. 2013). In the case of associating systems, constant ratio of signal/mg upon association is assumed but must be proven. The main concern with using labeled samples is that the label does not interfere with the molecular behavior. It is possible that label–label interactions occur in complexes that quench or enhance the signal intensity. This is not usually seen in tracer experiments; experiments with 100% labeled material are more likely to exhibit quenching behavior depending upon the locations or separation of the labels in oligomeric complexes. (An example of quenching in FDS data will be presented elsewhere.)

The nonideality equations used here (1–4) are referred to as phenomenological equations meaning empirical or experimental. Thermodynamic nonideality, BM1, refers to excluded volume (Tanford 1961) and charge effects as described by the Donnan term in osmotic pressure (Scatchard 1946) or sedimentation equilibrium experiments (Roark and Yphantis 1971). Hydrodynamic nonideality, ks, also refers to excluded volume and charge effects as reflected in the backflow during an SV experiment where displaced solvent must replace the volume a macromolecule vacates during sedimentation. The sedimenting particle has an effective solvated volume Vs that reflects hydration, shape through an f/fo term, and the entrained solvent captured within the Debye length (Fuoss 1959; Fuoss and Onsager 1932; Rowe 1977, 1992). It is worth noting that the effective charge is much smaller than typically assumed due to charge screening and anion binding (Laue 2011; Laue and Shire 2020). This in principle reflects the Stokes radius of the effective electroneutral sphere. We are aware that a good ks model system that controls for shape and charge would be useful here for investigating the dependence of ks on charge. Our purpose is not to assign theoretical values but rather to present the impact of typical empirical values for these parameters. It is more important that the impact of these values is graphically presented, and that their measurement by NLLS fitting to proper models, nonideality plus association, be clearly outlined.

It is generally assumed that the correct concentration scale for high concentration work is volume fraction (Ross and Minton 1977) and recent AUC studies on high concentration use the volume fraction Φ scale (Chaturvedi et al. 2018; Chaturvedi and Schuck 2019). The conversion is in principle a linear transformation, as simple as νc (Broide et al. 1991), but more realistically Vsc where the swollen or effective volume is used (Rowe 1992; Chaturvedi and Schuck 2019). At high volume fraction the shape and packing considerations may become important for mAbs (Garidel et al. 2017). The conversion of weight concentration to volume fraction can be verified from intrinsic viscosity measurements: \({\Phi }_{\mathrm{eff}}=[\eta ]c/2.5\), where 2.5 is for spheres (Cantor and Schimmel 1980). Therefore, we are currently making viscosity measurements on mAbs and serum proteins and will present those data in a subsequent publication along with AUC SV measurements. It has also been observed that mAbs deviate from linear viscosity behavior above approximately 50 mg/ml. (The Shire group has investigated the impact of association on viscosity and the problems it causes for drug delivery by injection (Liu et al. 2005; Yadav et al. 2011b, 2012).) We have confirmed these nonlinear viscosity observations and furthermore observe similar deviations in sw vs concentration plots. This has implications for SEDANAL analysis and requires a second order term, ks2c2, in the s/so and D/Do phenomenological Eqs. (3, 4). This and a higher order 3rd virial coefficient, CM1, were previously incorporated into SEDANAL. The polymer field refers to this nonlinear deviation as clustering or low energy attraction (~ 5 kT or 3 kcal/mol) near the boundary for liquid–liquid phase separation (Fiore et al. 2018). Clustering is clearly a distinct phenomena to the self-association we describe here.