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

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=ekxc then

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.

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.

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

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

where A>0 is the upper limit, kd−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)

depends only on the parameter d, and varies between 1e (limit value as d→1, the limit function is the Gompertz function) and 1 (limit value as d→∞), 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 1e.

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 ν=d−1instead of d of the Richards function when ν > 0 as in (13, 15, 25, 26):

Here A > 0 is the upper limit, kν>0 is the rate of the approximate exponential growth at minus infinity, and ν is a shape parameter determining the inflection rate value (ν+1)−1ν.

These functions are solutions of the differential equation.

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 (1+ν)1νA (13, 26):

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

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

when evaluating the measurement results. This is a concave function with an upper asymptote of lnA and a leftmost asymptote with a slope of kν.

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.

Its natural logarithm is

Both terms can be normalized to their own inflection points:

and

with the visual advantage of moving xi and yi to the same horizontal line.

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 ( lnFI) vs logarithm antigen concentration ( x) colored with the different serum dilutions, and the logarithm fluorescence intensity ( lnFI) vs serum dilution ( y) colored with the different antigen concentrations (Figure 4)

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

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.

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.

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.

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.

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.

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

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 (α).

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.

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.

where KD’ and KD are apparent and true equilibrium dissociation constants, [Ab]_s is the serum antibody concentration, ΔH° and Δ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.

The fluorescent signal intensity obtained by measuring bound [Ab], FI[Ab]bound, is then given by the Equation 12

where FIAb is the factor for converting bound [Ab] to fluorescence signals and can be obtained by calibration measurements.