Application of the Richards function to serum antibody titration

Conventional approaches to the titration of serum antibody binding use mid-point or end-point titers that are in a relative space and are therefore difficult to standardize. Here we propose the use of the thermodynamic titer, which, under appropriate measurement conditions, is a universal measure of the thermodynamic activity of serum antibodies. We show that the interpretation of the generalized logistic function as applied to biochemical binding events is possible using analogies to relative and absolute growth rates and size, which applies to the products of the reaction. Such deeper interpretation reveals the biological meaning of the asymmetry parameter of the function as a proportionality factor to ideal binding conditions. The use of a universally applicable and thermodynamically meaningful serum antibody titer could improve systematic mapping and understanding of antibody function.

1 Introduction

Antibody titration is similar to chemical titration in the sense that serially diluted serum is added to a reaction containing antigen. Unlike a chemical titration, the concentration of antibodies in serum is unknown and the concentration of antigen in the reaction is not necessarily known. As a consequence, the only way of expressing the result is a unitless index, the titer, which is expressed as the reciprocal of the dilution factor. The value of the titer for a serum antibody that reaches the predefined attributes when diluted 1:5000 is thus 5000. There are several approaches to characterize the titer (1). Endpoint titration considers the highest dilution where there is a positive signal (2). Midpoint (also called halfway) titration identifies the dilution factor where the signal falls to 50% of its starting value, often by fitting functions (3). By fitting a function to the binding curve one can interpolate to identify the dilution factor with greater precision. The use of antibody standards with known concentration allows transformation of titers into biochemical units of concentration. However, since monoclonal antibodies, their mixes or polyclonal antibody mixtures cannot precisely mimic a serum antibody binding curve these approaches are not accurate. The reason for this is the heterogeneity of an immune response: many different B-cell clones respond to an antigenic challenge and the composition of antibodies with respect to structure, affinity and concentration of these clones can be extremely variable.

The titration curve is a sigmoid growth curve that increases at an increasing rate until its inflection point and then increases at a decreasing rate approaching its upper limit. One of the first curves used to describe limited growth is the logistic curve, which is the graph of the logistic function first introduced by Verhulst (4). The various parametrizations of the family of logistic functions are known to immunologists as the four-parameter logistic function or 4PL (5, 6), as the Hill-equation for biochemists fitting ligand binding assays (7–9), and as the Langmuir-equation for modeling surface adsorption for biophysicists (10). From the theoretical point of view logistic functions model calibrated binary classification problems (11). Binary here refers to either being bound or unbound (free) and calibration means we assess the probability of these two states using a calibrated chemical force, the concentration of one of the reactants.

Because of the above discussed issues with standardization, most of the serological assays for specific antibodies express results in arbitrary units, which are not comparable across platforms, antigens and not even antibody isotypes. To circumvent the problem of obtaining results in incomparable arbitrary units we designed a method where antigen with known concentration is titrated simultaneously with serum antibody titration. This dual-titration approach built on the unique stoichiometry of microspot reactions and the practicality of running several reactions on the small surface area of a protein microarray (12). To fit the binding curves in two dimensions and to introduce flexibility into the binding curves of dual titration we applied the combination of two generalized logistic functions (13). Here we describe the rationale of our mathematical approach, show the analogies with growth curves and propose a general use for the Richards function in serum antibody titration.

2 Sigmoid functions and functions derived from sigmoid functions

2.1 Sigmoid functions

Sigmoid functions are increasing functions that have horizontal asymptotes at both ends and a single inflection point where the function is convex to the left and concave to the right. In applications, the left asymptote is usually zero. Many authors have summarized the sigmoid functions used in biological models and their properties (14–18, 25). Sigmoid functions are used, among others, to describe the size of finitely growing populations as a function of time (4) and the change in body sizes (19) or tumor sizes (20) also as a function of time. In immunochemistry, the concentration of the immune complex formed in antibody-antigen reactions is a sigmoid function as a function of both the logarithm of the antibody concentration and the logarithm of the antigen concentration (12, 21), in contrast to the exponential dependence obtained in the simple polyclonal model assuming constant antibody concentration (33) as a function of the logarithm of antibody concentration.

The definition of the sigmoid function only gives conditions for the shape of the graph of the function. We obtain more well-founded functions as solutions to the differential equation of the mathematical model that describes the phenomenon, because there is not only shape similarity between the experimental data curve and the graph of the function, but the graph also represents a solution to the differential equation that describes the dynamics of the phenomenon.

The most commonly used sigmoid functions in biology are the logistic function (22, 23), the Gompertz function (20, 24), and the Richards function (19). A special case of the Richards function is the logistic function, and its limiting case is the Gompertz function, and all of these are derived from differential equation models (13, 14, 16–18, 25). We mention that many continuous probability distribution functions are sigmoid, e.g. the normal and, for certain parameters, the Weibull.

The logistic function is

$\begin{array}{ll}L(x)=\frac{A}{1+b{e}^{-kx}}& (1)\end{array}$

where $A,b,k>0$, and its graph, the logistic curve, is symmetric about its inflection point. Here A is the limit at infinity or upper limit, k>0 is the approximate exponential growth rate at minus infinity, and b is the shift parameter. If we set $b=e^{k{x}_c}$ then

$$L(x)=\frac{A}{1+e^{-k(x-x_c)}}$$

and with this parametrization the real number x_c is the inflection point of the function. This family of functions is the general solution of the logistic differential equation.

$$\frac{dL}{dx}=kL(1-\frac{L}{A})$$

with a range of values between 0 and A. If an additive constant is included, we get the four-parameter logistic function or 4PL as a function of e^x which allows for a non-zero left asymptote (see the Supplementary Materials).

We get a natural generalization of the logistic function (Equation 1) by raising its denominator to a power, which is called the Richards function or generalized logistic function.

$$R(x)=\frac{A}{{(1+b{e}^{-kx})}^m}$$

For this function, $A>0$ is the upper limit, and for $k,m>0$, km is the rate of the initial approximately exponential growth. Unfortunately, the function given in the form

$$R(x)=\frac{A}{{(1+e^{-k(x-x_c)})}^m}$$

has an inflection point at x_c only if $m=1$. Therefore, a commonly used form of the Richards function for $m>0$ is

$\begin{array}{ll}R(x)=A{(1+(d-1)e^{-k(x-x_i)})}^{\frac{1}{1-d}}& (2)\end{array}$

where A>0 is the upper limit, $\frac{k}{d-1}>0$ is the rate of the approximate exponential growth at minus infinity, d>1 is the asymmetry parameter and the real number x_i is the inflection point of the function (Figure 1). The ratio of the value at the inflection point and the upper limit (in short inflection rate value)

Figure 1

Figure 1. Richards growth function on the lin–lin scale (A) and the log–lin scale (B) for different shape parameters d, with all other parameters held constant. The x-axis represents the independent variable x (e.g., logarithm antigen concentration or logarithm serum dilution). The dots denote the inflection points for each curve.

$$\frac{R(x_i)}{A}=d^{\frac{1}{1-d}}$$

depends only on the parameter d, and varies between $\frac{1}{e}$ (limit value as $d\to 1$, the limit function is the Gompertz function) and 1 (limit value as $d\to \infty$), and d = 2 returns the logistic function.

For 0 < d < 1, the Richards function (Equation 2) is usually defined starting from its zero point, in an interval infinite from the right. Its inflection rate value changes between 0 and $\frac{1}{e}$.

If we add an arbitrary constant to the Richards function, we obtain the five-parameter logistic function 5PL as a function of e^x (see the Supplementary Materials).

In our work, we use the parameter $\nu =d-1$instead of d of the Richards function when ν > 0 as in (13, 15, 25, 26):

$\begin{array}{ll}R(x)=A{(1+\nu e^{-k(x-x_i)})}^{-\frac{1}{\nu }}& (3)\end{array}$

Here A > 0 is the upper limit, $\frac{k}{\nu }>0$ is the rate of the approximate exponential growth at minus infinity, and ν is a shape parameter determining the inflection rate value ${(\nu +1)}^{-\frac{1}{\nu }}$.

These functions are solutions of the differential equation.

$$\frac{dR}{dx}=\frac{k}{\nu }R(1-{(\frac{R}{A})}^{\nu })$$

The slope, or rate of growth, of a sigmoid function is greatest at its inflection point. Since this point is an important location also for the Richards function, it is reasonable to normalize the function at the inflection point, i.e., to achieve a value of 1 at the inflection point using the appropriate constant multiplier $\frac{{(1+\nu )}^{\frac{1}{\nu }}}{A}$ (13, 26):

$$R_{\nu }(x)={(\frac{1+\nu }{1+\nu e^{-k(x-x_i)}})}^{\frac{1}{\nu }}$$

For various parametrizations and other properties of the sigmoid functions discussed above, see the Supplementary Materials.

2.2 Functions derived from sigmoid functions

In our work, we take the logarithmic version of the Richards function, Equation 3, which yields

$\begin{array}{ll}\text{ln}R(x)=\text{ln}A-\frac{1}{\nu }\text{ln}(1+\nu e^{-k(x-x_i)})& (4)\end{array}$

when evaluating the measurement results. This is a concave function with an upper asymptote of lnA and a leftmost asymptote with a slope of $\frac{k}{\nu }$.

A double Richards function is a two-variable function that is the product of a Richards function of each of the two variables. A common rate parameter k is used for both terms to represent identical slope behavior, consistent with the symmetry assumed in our experimental framework.

$\begin{array}{ll}R(x,y)=R_1(x)R_2(y)=C{(1+{\nu }_1{e}^{-k(x-x_i)})}^{- \frac{1}{{\nu }_1}}{(1+{\nu }_2{e}^{-k(y-y_i)})}^{- \frac{1}{{\nu }_2}}& (5)\end{array}$

Its natural logarithm is

$$\ln R(x,y)=\text{ln}(R_1(x)R_2(y))=\text{ln}C-\frac{1}{{\nu }_1}\text{ln}(1+{\nu }_1{e}^{-k(x-x_i)})-\frac{1}{{\nu }_2}\text{ln}(1+{\nu }_2{e}^{-k(y-y_i)})$$

Both terms can be normalized to their own inflection points:

$\begin{array}{ll}{R}_n(x,y)=R_1(x)R_2(y)=C_n{(\frac{1+{\nu }_1}{1+{\nu }_1{e}^{-k(x-x_i)}})}^{\frac{1}{{\nu }_1}}{(\frac{1+{\nu }_2}{1+{\nu }_2{e}^{-k(y-y_i)}})}^{\frac{1}{{\nu }_2}}& (6)\end{array}$

and

$\begin{array}{ll}\ln R_n(x,y)=\text{ln}(R_1(x)R_2(y))=\text{ln}{C}_n+\frac{1}{{\nu }_1}\text{ln}(\frac{1+{\nu }_1}{1+{\nu }_1{e}^{-k(x-x_i)}})+\frac{1}{{\nu }_2}\text{ln}(\frac{1+{\nu }_2}{1+{\nu }_2{e}^{-k(y-y_i)}})& (7)\end{array}$

with the visual advantage of moving $x_i$ and $y_i$ to the same horizontal line.

3 Graphical presentation of Richards curves

Figure 2 shows how the Richards function (Equation 3) behaves as its parameters ($A$, $\nu$, $x_i$ and $k$) are varied individually. Increasing $A$ proportionally shifts the function vertically, raising its upper asymptote. The asymmetry parameter ν controls the steepness of the curve: higher values of $\nu$ give less steep transitions between the asymptotes. Varying the inflection point $x_i$ results in a horizontal shift of the curve, i.e., increasing $x_i$ moves the entire sigmoid curve to the right. When $A$ and $\nu$ remain the same, then regardless of the value of $k$, $R(x_i)$ will be the same. A larger $k$ implies a more rapid change, it determines how gradual the transition of the sigmoid curve is. Similar presentation is seen in (5).

Figure 2

Figure 2. The Richards function (Equation 3) plotted as each of its parameters – $A$ (A), $\nu$ (B), $x_i$ (C) and $k$ (D) – is varied individually. The fixed parameter values are $A=8000$, $\nu =0.5$, $x_i=-13$ and $k=1$, while the varying parameter values are specified in each corresponding subplot. The x-axis range (–15 to –11) corresponds to the range of log transformed concentrations used in the experimental dataset introduced later in Section 3.

Similarly Figure 3 shows how the Richards function (Equation 3) behaves on the logarithmic scale (Equation 4) as its parameters ($lnA$, $\nu$, $x_i$ and $k$) are varied individually. Increasing $lnA$ shifts the function vertically, raising its upper asymptote. The asymmetry parameter ν controls the steepness of the curve: higher values of $\nu$ produce curves that approach the asymptotes more slowly on one side. Varying the inflection point $x_i$ results in a horizontal shift of the curve, while changes in $k$ determines how gradual the transition of the sigmoid curve is.

Figure 3

Figure 3. The Richards function (Equation 3) plotted on the logarithmic scale (Equation 4) as each of its parameters – $lnA$ (A), $\nu$ (B), $x_i$ (C) and $k$ (D) – varied individually. The fixed parameter values are $lnA=9$, $\nu =0.5$, $ln{x}_i=-13$ and $k=1$. The varying parameter values are specified in each corresponding subplot. The x-axis range (–15 to –11) corresponds to the range of log transformed concentrations used in the experimental dataset introduced later in Section 3.

3.1 Illustration using experimental data

To illustrate the Richards model applicability we fitted it to experimental data obtained from a representative serum sample PS332, which was part of the commercially available set of serum samples from individuals with PCR confirmed SARS-CoV-2 infection from Raybiotech (Peach Tree corners, Georgia, USA). IgG reactivity against the RBD domain of the virus was measured as described previously (27). The reactivity of IgG against the tested viral antigen was confirmed using the experimental technology that builds on our model and simultaneously fits two Richards functions to the measurement data. The example is intended for methodological demonstration rather than biological generalization.

Data are best visualized by plotting the logarithm of fluorescence intensity ($\ln FI$) vs logarithm antigen concentration ($x$) colored with the different serum dilutions, and the logarithm fluorescence intensity ($\ln FI$) vs serum dilution ($y$) colored with the different antigen concentrations (Figure 4)

Figure 4

Figure 4. Experimental data (IgG 332) is presented in two dimensions. Logarithm of fluorescence intensity ($\ln FI$) vs logarithm antigen concentration ($x$) colored with different serum dilutions (A) and logarithm of fluorescence intensity ($\ln FI$) vs logarithm serum dilution ($y$) colored with different antigen concentrations (B).

In addition, these data can also be presented in 3 dimensions as shown in Figure 5.

Figure 5

Figure 5. Experimental data (IgG 332) is presented in three dimensions. Measurements of logarithm of fluorescence intensity ($\ln FI$) vs logarithm antigen concentration and logarithm serum dilution colored with different serum dilutions (A) and colored with different antigen concentrations (B).

3.2 Statistical performance

The estimation of the Richards parameters is feasible with standard nonlinear least squares methods, and it generally provides robust fits especially when the full S-shaped curve is well sampled. However, convergence issues can arise when data are sparse or when the inflection point is poorly sampled, which can affect the stability and reliability of parameter estimates. Since the literature on the statistical properties of the single Richards curve is already extensive, we now focus on the statistical performance of the double Richards model.

3.3 Statistical performance of the double Richards model

Unlike simpler sigmoidal models, the double Richards model (Equation 5) requires the estimation of two sets of nonlinear parameters, which presents both computational and interpretive challenges. We applied nonlinear regression to fit the double Richards model and additionally, we incorporated Bayesian inference to improve the estimating procedure.

Fitting the double Richards model relies on iterative nonlinear optimization algorithms, typically the Levenberg-Marquardt method, which poses significant statistical challenges. Achieving convergence in nonlinear regression is not guaranteed, particularly for complex models like the double Richards, where parameter spaces are large and often have flat or multimodal likelihood surfaces.

In our experience, naive initialization of parameters frequently led to convergence failure or biologically implausible estimates. To improve convergence reliability, we implemented a data-driven procedure to generate starting values tailored to the structure of the double Richards function. Central parameters were set using medians or relevant quantiles, depending on the skewness of the variables, while shape parameters were set to 1. This approach provided reasonable starting values across our datasets and helped mitigate common convergence failures.

We also considered other techniques such as using derivatives and piecewise regression but even with these procedures convergence was not universal, especially in datasets with high noise levels or limited range. However, we also note that having good initial parameter estimates does not guarantee convergence in nonlinear regression. Even when starting values are close to the true parameters, convergence may fail due to the complex (often flat) shape of the likelihood surface.

Despite obtaining good starting values, convergence was not guaranteed due to the model’s complexity. Therefore, we also employed a Bayesian estimation framework using weakly informative priors. Parameter inference was performed via Markov Chain Monte Carlo (MCMC) sampling, which provided greater stability and more reliable estimates, especially when the likelihood surface was complex or multi-modal.

3.4 Illustration of the double Richards model using experimental data

To demonstrate the applicability and performance of the double Richards model, we fitted it to the experimental biding data (serum sample PS332 from an individual with PCR confirmed SARS-CoV-2 infection, having measured the IgG reactivity against the RBD domain of the virus) introduced in Section 3.1. Parameter estimates were obtained using both nonlinear least squares and Bayesian methods as shown in Table 1.

Table 1

Table 1. Parameter estimation from the double Richards model (Equation 7) are given using both nonlinear least squares and Bayesian methods using the IgG PS332 data.

To demonstrate the performance of the double Richards model in fitting complex sigmoidal curves, we applied it to experimental binding data obtained from a serum sample (PS332). This sample was derived from an individual with PCR-confirmed SARS-CoV-2 infection, and the IgG antibody reactivity was measured against the receptor-binding domain (RBD) of the virus. The data, introduced in Section 3.1, displays a biphasic response characteristic, making it a suitable candidate for modeling with a double-sigmoid function.

Figure 6 shows the estimated fluorescence intensity based on the above model.

Figure 6

Figure 6. The double Richards function (Equation 6) is plotted on the logarithm scale (Equation 7) with parameters listed in Table 1, shown together with the corresponding experimental data (IgG 332). Logarithm of fluorescence intensity (ln FI) vs logarithm of antigen concentration (x) colored by serum dilution (A) and logarithm of fluorescence intensity (ln FI) vs logarithm of serum dilution (y) colored by antigen concentration (B).

To further illustrate the behavior of the fitted model, these data are also presented in 3-dimensions (Figure 7), showing how the modeled response (logarithm of fluorescence intensity) changes concurrently with both independent variables - the logarithm of antigen concentration and the logarithm of serum dilution. This representation highlights the interaction between the two variables and the overall shape of the fitted surface.

Figure 7

Figure 7. Three-dimensional plot of the double Richards function (Equation 6) on a logarithmic scale (Equation 7), showing the logarithm of fluorescence intensity as a function of the logarithm of antigen concentration and logarithm of serum dilution. Parameter values were estimated from the experimental data (Table 1).

Experimental data (IgG 332) is presented in three dimensions. Measurements of logarithm of fluorescence intensity ($\ln FI$) vs logarithm antigen concentration and logarithm serum dilution colored with different serum dilutions (A) and colored with different antigen concentrations (B).

4 Application of logistic functions to biochemical reactions

4.1 Analogy to growth curves

Sigmoid growth curves describe cumulative processes: cells accumulate forming increasing mass, trees accumulate creating forests, individuals accumulate resulting in growth of populations (14, 15, 28, 29). Such biological growth is observed as a function of time. During biochemical reactions products of the reaction can also accumulate, their concentration grows, however, titration curves are expressed as a function of chemical forces instead of time. In our case, when the biochemical reaction is the binding of Ab and Ag molecules, these chemical forces are expressed as the logarithms of the concentration of the titrated component. The reason for taking the logarithm is to obtain quantities related to chemical energy or chemical potential. The observation of growth here therefore does not mean the observing of change of product concentrations of a given reaction in time but rather the equilibrium product concentrations from several distinct reactions where different chemical energies were applied. In each of these different measurements the reaction has to reach equilibrium in order to make them comparable and suitable for analysis. This “growth curve” is the binding curve, with its discrete measurement points in the dimension of chemical potential instead of time.

When modelling growth with logistic functions, two important and useful transformations of the growth function can be made (30). The derivative of the function, which describes absolute growth rate, tells us how fast the concentration of products is increasing as we increase the chemical force applied. Even more useful is the relative growth rate, which is described by the derivative of the logarithm of the growth function. Relative growth rate is the rate of change of size compared to current size. This is a monotonously decreasing function, as the size increases the relative rate of growth slows down. This is because the resources supporting the growing entity gradually become smaller and smaller in comparison to the needs of the entity. For a binding reaction that is interpreted as the waning of chemical energy supporting the reaction against the energies counteracting product formation. In other words, this function represents how the titrated chemical energy is used up for supporting the reaction (Figure 8). It is a mathematical identity that the product of the growth function and the relative growth rate function yields the absolute growth rate function. In thermodynamic terms that corresponds to the product of concentration ([Ab]bound) and activity coefficient (γ), yielding the relative thermodynamic activity (α).

Figure 8

Figure 8. Relationship between curves of growth (A-C) and biochemical reaction (D-F). Analogies are apparent between pairs (A, D, B, E, C, F) of related phenomena.

4.2 Interpretation of deformation of the logistic curve

Growth curves are not always as ideal as would be suggested by a logistic growth curve (Figure 9A). The generalization of the logistic function allows for deformations of the curve by the introduction of another parameter into the function: the asymmetry parameter. Non-ideality of binding curves has also been observed for immunoassays, where the fitting of the curves was also improved by adding another parameter to the function, transforming four-parameter logistic function (4PL) into five-parameter logistic function (5PL) (5, 31, 32). The introduction of the novel parameter yields functions that are solutions of the Richards differential equation, but the parametrization of these solutions has important consequences for the interpretation of biochemical meaning of the parameters (see (13, 33)).

From the molecular point of view the deformation of the binding curve means that the accumulation of products (proportional to the absorbed chemical energy) and the decrease in available energy are not symmetric: we are “losing” energy somewhere. This loss is only in the observation of course, and suggests that our model does not count with factors that contribute to the binding. Such an asymmetry in the measurement can arise when the measurement system itself is biased for the detection of certain forms of energy. Microspot immunoassays are also known as mass independent immunoassays because under appropriate conditions the total number of molecules (mass) is not affecting the measurement, only the concentration (thermodynamic activity). The advantages of mass independence can be exploited for antibody sandwich type measurements of antigen analytes (34, 35) and also for antigen specific serum antibody measurements (12). In a stricter thermodynamic sense microspot assays can allow the separation of the measurement of two thermodynamic potentials: enthalpy and entropy (Figures 9B, C) (36, 37). By varying the density of surface immobilized antigen molecules, we can adjust the number of binding sites for non-covalent bonding per unit area. That way we control the enthalpic forces contributing to binding. On the other hand, by varying the concentration of antibody molecules in solution, we adjust the number of ‘particles’ contributing to the filling of energy levels of the system. By particles here we mean antibody molecules, which populate the binding energy landscape, and which attain different energy levels by conformational changes. Thus, when titrating surface antigen density or antibody concentration in solution, we decouple the measurement of contribution of enthalpic and entropic forces to the binding.

Figure 9

Figure 9. Deformation of ideal curves by observation bias. When product formation is a function of free energy (A), then the titration of enthalpic (B) or entropic (C) components of binding reaction will result in deformations of the logistic curve. While Ag concentration can be varied over the relevant reasonable range (B), serum titration starts with undiluted serum, represented by filled circles as end-points (C).

4.3 Two Richards functions for two thermodynamic potentials

By decoupling the measurement of enthalpic and entropic contributions to binding, we can estimate the partial contribution of these thermodynamic potentials. Because these two chemical forces jointly contribute to binding, the estimation of the partial contribution of one results in solving the partial contribution of the other, as well. We can incorporate the deformation caused by enthalpy or entropy-biased measurements (Figure 10) into the asymmetry parameter of the Richards function as defined by Equations 8, 9.

Figure 10

Figure 10. Effects of asymmetry parameters on normalized functions. The growth of bound Ab concentration, as a function of reaction enthalpy (A) and entropy (B). Apparent affinity under these conditions is always greater than true affinity, (KD’< KD), therefore ln(${\nu }_{Ag}$) is always negative. Excess entropy can both hinder or enhance binding, ln(${\nu }_{Ab}$) can therefore be both positive or negative, corresponding to the observed total serum [Ab] being higher or lower than KD.

$\begin{array}{ll}\text{ln}({\nu }_{Ag})=\text{ln}(\frac{KD\text{'}}{KD})=(\Delta H-\Delta H°)/RT=\Delta \Delta H/RT& (8)\end{array}$

$\begin{array}{ll}\text{ln}({\nu }_{Ab})=\text{ln}(\frac{{[Ab]}_s}{KD})=(\Delta S-\Delta S°)/R=\Delta \Delta S/R& (9)\end{array}$

where KD’ and KD are apparent and true equilibrium dissociation constants, [Ab]_s is the serum antibody concentration, $\Delta H°$ and $\Delta S°$ stand for ideal, equal contribution of enthalpy and entropy to binding, ΔH and ΔS are total enthalpic and entropic contributions to binding, while ΔΔH and ΔΔS are enthalpic and entropic parts of excess Gibbs free energy, respectively (38). Energies are expressed with reference to the universal gas constant R and thermodynamic temperature T.

Whilst the two asymmetry parameters could also be estimated separately, technically we can simultaneously measure and estimate the two contributions. This is carried out by creating a measurement matrix, where rows and columns correspond to binding signals obtained by titrating enthalpic and entropic forces, namely ln[Ag] and ln[Ab]. It is due to the measurement conditions that the titration of Ag and Ab reflect distinct chemical forces: Ag molecules are immobilized on a surface in the presence of huge excess of Ab in solution (12, 27). We can then fit two Richards curves with asymmetry parameters, defined by Equations 10, 11, simultaneously, one for Ag and one for Ab titration.

$\begin{array}{ll}{R}_{1,\nu }([Ag])~{(\frac{1+{\nu }_{Ag}}{1+{\nu }_{Ag}{e}^{-k(\text{ln}[Ag]-\text{ln}[Ag]°)}})}^{\frac{1}{{\nu }_{Ag}}}& (10)\end{array}$

$\begin{array}{ll}{R}_{2,\nu }([Ab])~{(\frac{1+{\nu }_{Ab}}{1+{\nu }_{Ab}{e}^{-k^{\prime }(\text{ln}[Ab]-\text{ln}[Ab]°)}})}^{\frac{1}{{\nu }_{Ab}}}& (11)\end{array}$

The fluorescent signal intensity obtained by measuring bound [Ab], $F{I}_{{[Ab]}_{bound}}$, is then given by the Equation 12

$\begin{array}{ll}F{I}_{{[Ab]}_{bound}}=F{I}_{Ab}*{[Ab]}_{bound}^{°}*R_{1,\nu }([Ag])*R_{2,\nu }([Ab])& (12)\end{array}$

where $F{I}_{Ab}$ is the factor for converting bound [Ab] to fluorescence signals and can be obtained by calibration measurements.

5 Practical guide for using Richards curves for titration

Most of the immune assays in use today employ only a single Ag concentration. Therefore, if the condition of massive excess of Ab in solution is met by the assay, it is still possible to fit an Ab titration curve without an Ag titration curve. While the affinity cannot be addressed without Ag titration it is possible to obtain universally comparable values characterizing serum antibodies: the dilution factor of the undiluted serum from the point of inflection, given by ${\nu }_{Ab}$ is the thermodynamic titer, also called thermodynamic concentration, [Ab]s/K_D. The thermodynamic titer is an indicator of the ability of antibodies to saturate Ag. When its value is 1, Ag is half saturated. Above one saturation increases, below one it decreases below 50%. This value is therefore independent of the affinity of serum antibodies.

By calibrating the signals representing Ab binding we can convert signal intensity to bound Ab concentration, [Ab]_b. This approach is used in one-point assays for obtaining results with concentration units, but it should be noted that the proper unit in that case is [Ab]_eq, referring to the fact that the observed concentration is equivalent to the Ab (often a monoclonal Ab) used for calibrating binding signals. We used known concentrations of serum Ab classes printed on microarrays as calibration curves (12, 27). Thus, while the thermodynamic concentration can be obtained without calibration, determination of the standard binding site occupancy, ${[Ab]}_{bound}^{°}$, requires calibration (Figure 11).

Figure 11

Figure 11. Comparison of one-point, titration and dual-titration assay outputs. Standard binding site occupancy ${[Ab]}_{bound}^{°}$, K_D, equilibrium dissociation constant.

The added advantage of carrying out a titration for two assay components, serum [Ab] and assay [Ag] simultaneously is that potential biases in the measurement of one can be compensated by the other. The combination of two Richards curves into a single fitting algorithm will result in a search for the intersection of the two curves at their respective inflection points, adjusting the three-dimensional surface defined by the combination of the two functions. As this is the point where both functions are changing at the fastest rate, a sensitive fitting can be obtained. Therefore, we suggest the use of dual-titration over the titration ranges expected to contain their inflection points.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.

Author contributions

ÁK: Writing – review & editing, Conceptualization, Writing – original draft, Visualization. KP: Writing – review & editing, Visualization. JP: Visualization, Conceptualization, Writing – review & editing, Writing – original draft. TP: Writing – review & editing, Conceptualization, Writing – original draft.

Funding

The author(s) declare financial support was received for the research and/or publication of this article. ÁK and TP was funded by the National Research, Development and Innovation Office in Hungary (RRF-2.3.1-21-2022-00006).

Conflict of interest

Authors KP and JP are employed by company Diagnosticum Zrt.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1657633/full#supplementary-material

References

1. Farnsworth CW, Case JB, Hock K, Chen RE, O’Halloran JA, Presti R, et al. Assessment of serological assays for identifying high titer convalescent plasma. Transfusion. (2021) 61(9):2658–2667. doi: 10.1111/trf.16580

PubMed Abstract | Crossref Full Text | Google Scholar

2. Frey A, Di Canzio J, and Zurakowski D. A statistically defined endpoint titer determination method for immunoassays. J Immunol Methods. (1998) 221:35–41. doi: 10.1016/S0022-1759(98)00170-7

PubMed Abstract | Crossref Full Text | Google Scholar

3. Phiri K and Grill L. Development of a candidate TMV epitope display vaccine against SARS-coV-2. Vaccines (Basel). (2024) 12. doi: 10.3390/vaccines12050448

PubMed Abstract | Crossref Full Text | Google Scholar

4. Verhulst P-F. Notice sur la loi que la population suit dans son accroissement. Correspondance mathématique physique. (1838) 10:113–21.

Google Scholar

5. Gottschalk PG and Dunn JR. The five-parameter logistic: a characterization and comparison with the four-parameter logistic. Anal Biochem. (2005) 343:54–65. doi: 10.1016/j.ab.2005.04.035

PubMed Abstract | Crossref Full Text | Google Scholar

6. Cumberland WN, Fong Y, Yu X, Defawe O, Frahm N, and De Rosa S. Nonlinear calibration model choice between the four and five-parameter logistic models. J Biopharm Stat. (2015) 25:972–83. doi: 10.1080/10543406.2014.920345

PubMed Abstract | Crossref Full Text | Google Scholar

7. Goutelle S, Maurin M, Rougier F, Barbaut X, Bourguignon L, Ducher M, et al. The Hill equation: a review of its capabilities in pharmacological modelling. Fundam Clin Pharmacol. (2008) 22:633–48. doi: 10.1111/j.1472-8206.2008.00633.x

PubMed Abstract | Crossref Full Text | Google Scholar

8. Buchwald P. A single unified model for fitting simple to complex receptor response data. Sci Rep. (2020) 10:13386. doi: 10.1038/s41598-020-70220-w

PubMed Abstract | Crossref Full Text | Google Scholar

9. Gesztelyi R, Zsuga J, Kemeny-Beke A, Varga B, Juhasz B, and Tosaki A. The Hill equation and the origin of quantitative pharmacology. Arch Hist Exact Sci. (2012) 66:427–38. doi: 10.1007/s00407-012-0098-5

Crossref Full Text | Google Scholar

10. Islam MA, Chowdhury MA, Mozumder MSI, and Uddin MT. Langmuir adsorption kinetics in liquid media: interface reaction model. ACS Omega. (2021) 6:14481–92. doi: 10.1021/acsomega.1c01449

PubMed Abstract | Crossref Full Text | Google Scholar

11. Kim S-C, Arun AS, Ahsen ME, Vogel R, and Stolovitzky G. The Fermi-Dirac distribution provides a calibrated probabilistic output for binary classifiers. Proc Natl Acad Sci USA. (2021) 118. doi: 10.1073/pnas.2100761118

PubMed Abstract | Crossref Full Text | Google Scholar

12. Papp K, Kovács Á, Orosz A, Hérincs Z, Randek J, Liliom K, et al. Absolute quantitation of serum antibody reactivity using the Richards growth model for antigen microspot titration. Sensors. (2022) 22. doi: 10.3390/s22103962

PubMed Abstract | Crossref Full Text | Google Scholar

13. Prechl J, Papp K, Kovács Á, and Pfeil T. The binding landscape of serum antibodies: how physical and mathematical concepts can advance systems immunology. Antibodies (Basel). (2022) 11. doi: 10.3390/antib11030043

PubMed Abstract | Crossref Full Text | Google Scholar

14. Tsoularis A and Wallace J. Analysis of logistic growth models. Math Biosci. (2002) 179:21–55. doi: 10.1016/S0025-5564(02)00096-2

PubMed Abstract | Crossref Full Text | Google Scholar

15. Tjørve E and Tjørve KMC. A unified approach to the Richards-model family for use in growth analyses: why we need only two model forms. J Theor Biol. (2010) 267:417–25. doi: 10.1016/j.jtbi.2010.09.008

PubMed Abstract | Crossref Full Text | Google Scholar

16. Tjørve KMC and Tjørve E. A proposed family of Unified models for sigmoidal growth. Ecol Modell. (2017) 359:117–27. doi: 10.1016/j.ecolmodel.2017.05.008

Crossref Full Text | Google Scholar

17. France J, Dijkstra J, and Dhanoa MS. Growth functions and their application in animal science. Anim Res. (1996) 45:165–74. doi: 10.1051/animres:19960637

Crossref Full Text | Google Scholar

18. Marusic M. Mathematical models of tumor growth. Zagreb, Croatia: Zagreb, Croatia: Department of Mathematics, University of Zagreb, Bijeniˇcka (1996)30.

Google Scholar

19. Richards FJ. A flexible growth function for empirical use. J Exp Bot. (1959) 10:290–301. doi: 10.1093/jxb/10.2.290

Crossref Full Text | Google Scholar

20. Gompertz B. On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of Life Contingencies. Philos Trans R Soc London. (1825) 115:513–83.

Google Scholar

21. Beatty JD, Beatty BG, and Vlahos WG. Measurement of monoclonal antibody affinity by non-competitive enzyme immunoassay. J Immunol Methods. (1987) 100:173–9. doi: 10.1016/0022-1759(87)90187-6

PubMed Abstract | Crossref Full Text | Google Scholar

22. Kuriyan J, Konforti B, and Wemmer D. The molecules of life: physical and chemical principles. Garland Sci. (2012). doi: 10.1201/9780429258787

Crossref Full Text | Google Scholar

23. Hingley PJ and Ouldridge EJ. The use of a logistic model for the quantitative interpretation of indirect sandwich enzyme labelled immunosorbent assays (ELISA) for antibodies and antigens in foot and mouth disease. Comput Biol Med. (1985) 15:137–52. doi: 10.1016/0010-4825(85)90027-7

PubMed Abstract | Crossref Full Text | Google Scholar

24. Vaghi C, Rodallec A, Fanciullino R, Ciccolini J, Mochel JP, Mastri M, et al. Population modeling of tumor growth curves and the reduced Gompertz model improve prediction of the age of experimental tumors. PloS Comput Biol. (2020) 16:e1007178. doi: 10.1371/journal.pcbi.1007178

PubMed Abstract | Crossref Full Text | Google Scholar

25. Ukalska J and Jastrzębowski S. Sigmoid growth curves, a new approach to study the dynamics of the epicotyl emergence of oak. Folia Forestalia Polonica. (2019) 61:30–41. doi: 10.2478/ffp-2019-0003

Crossref Full Text | Google Scholar

26. Beg C, Byeon J, Berman N, and Kieffer J. Correlative activation free energy (CAFE) model: Application to viscous processes in inorganic oxide glass-formers. Acta Mater. (2025) 283:120538. doi: 10.1016/j.actamat.2024.120538

Crossref Full Text | Google Scholar

27. Kovács Á, Hérincs Z, Papp K, Kaczmarek JZ, Larsen DN, Stage P, et al. In-depth immunochemical characterization of the serum antibody response using a dual-titration microspot assay. Front Immunol. (2025) 16:1494624. doi: 10.3389/fimmu.2025.1494624

PubMed Abstract | Crossref Full Text | Google Scholar

28. Birch C. A new generalized logistic sigmoid growth equation compared with the Richards growth equation. Ann Bot. (1999) 83:713–23. doi: 10.1006/anbo.1999.0877

Crossref Full Text | Google Scholar

29. Yin X, Goudriaan J, Lantinga EA, Vos J, and Spiertz HJ. A flexible sigmoid function of determinate growth. Ann Bot. (2003) 91:361–71. doi: 10.1093/aob/mcg029

PubMed Abstract | Crossref Full Text | Google Scholar

30. Pommerening A and Muszta A. Relative plant growth revisited: Towards a mathematical standardisation of separate approaches. Ecol Modell. (2016) 320:383–92. doi: 10.1016/j.ecolmodel.2015.10.015

Crossref Full Text | Google Scholar

31. Dawson DA, Genco N, Bensinger HM, Guinn D, Il’giovine ZJ, Wayne Schultz T, et al. Evaluation of an asymmetry parameter for curve-fitting in single-chemical and mixture toxicity assessment. Toxicology. (2012) 292:156–61. doi: 10.1016/j.tox.2011.12.006

PubMed Abstract | Crossref Full Text | Google Scholar

32. Findlay JWA and Dillard RF. Appropriate calibration curve fitting in ligand binding assays. AAPS J. (2007) 9:E260–7. doi: 10.1208/aapsj0902029

PubMed Abstract | Crossref Full Text | Google Scholar

33. Pfeil T and Herbály B. A linear model for polyclonal antibody–antigen reactions. Math Comput Simul. (2022) 198:20–30. doi: 10.1016/j.matcom.2022.02.004

Crossref Full Text | Google Scholar

34. Ekins RP. Multi-analyte immunoassay. J Pharm BioMed Anal. (1989) 7:155–68. doi: 10.1016/0731-7085(89)80079-2

PubMed Abstract | Crossref Full Text | Google Scholar

35. Saviranta P, Okon R, Brinker A, Warashina M, Eppinger J, and Geierstanger BH. Evaluating sandwich immunoassays in microarray format in terms of the ambient analyte regime. Clin Chem. (2004) 50:1907–20. doi: 10.1373/clinchem.2004.037929

PubMed Abstract | Crossref Full Text | Google Scholar

36. Chu X, Gan L, Wang E, and Wang J. Quantifying the topography of the intrinsic energy landscape of flexible biomolecular recognition. Proc Natl Acad Sci USA. (2013) 110:E2342–51. doi: 10.1073/pnas.1220699110

PubMed Abstract | Crossref Full Text | Google Scholar

37. Liu S-Q, Ji X-L, Tao Y, Tan D-Y, Zhang K-Q, and Fu Y-X. Protein folding, binding and energy landscape: A synthesis. In: Kaumaya P, editor. Protein engineering. InTech. (2012) doi: 10.5772/30440

Crossref Full Text | Google Scholar

38. Prechl J. Statistical thermodynamics of self-organization of the binding energy super-landscape in the adaptive immune system. arXiv. (2024).

Google Scholar