Sicegar: R package for sigmoidal and double-sigmoidal curve fitting

Sigmoidal and double-sigmoidal models

The simplest sigmoidal curve can be uniquely determined by three parameters, the maximum value, the midpoint, and the slope of the curve (Fig. 1A). To represent such a sigmoidal curve, the sicegar package uses the logistic function (Verhulst, 1845) $$I(t)=f_{\text{sig}}(t)=\frac{I_{\text{max}}}{1+\text{exp}(-a_1(t-t_{\text{mid}}))}.$$ (1) Here, I(t) is the intensity time course, given as a function of time t. The three parameters to be fitted are I_max, t_mid, and a₁. The parameter I_max represents the maximum intensity observed, the parameter t_mid indicates the time at which intensity has reached half of its maximum, and the parameter a₁ is related to the slope of I(t) at t = t_mid via the formula d/dt I(t)t = t_mid = a₁I_max/4 (Fig. 1A).

A double-sigmoidal curve will generally need at least six parameters, two midpoints and slopes parameters each, a maximum value, and an asymptotic final value after decay (Fig. 1B). We define the double-sigmoidal model by combining two regular sigmoidal functions. This combined function rises from 0 to a maximum value I_max at time t^* and then decays to the final value I_final. The function is divided into two parts, one to the left of the maximum, for t < t^* (the growth phase) and one to the right of the maximum, for t > t^* (the decay phase). The two parts are rescaled separately, such that the maximum value of I(t) is given by I_max and the final value of I(t) for large times is given by I_final > 0.

To define our double-sigmoidal model, we begin with a simple double-sigmoidal function, obtained by multiplying together two sigmoidal functions: $$f_{\text{dsig-base}}(t)=\frac{1}{1+\exp (-{a^{\prime }}_1(t-{t^{\prime }}_{\text{mid1}}))}\frac{1}{1+\exp (-{a^{\prime }}_2(t-{t^{\prime }}_{\text{mid2}}))}.$$ (2)

We then define the time t^* at which this function is maximal, $$t^{*}=\arg \max f_{\text{dsig-base}}(t),$$ (3) and the value of the function at that time point, $$f_{\text{max}}=f_{\text{dsig-base}}(t^{*}).$$ (4)

We then write $$I(t)=f_{\text{dsig}}(t)=\{\begin{array}{ll}{c}_1{f}_{\text{dsig-base}}(t)\hfill & \text{for}t\le t^{\text{*}}(\text{growth}\text{phase})\hfill \\ c_2{f}_{\text{dsig-base}}(t)+I_{\text{final}}\hfill & \text{for}t>t^{\text{*}}(\text{decay}\text{phase})\hfill \end{array},$$ (5) with $$c_1=\frac{I_{\text{max}}}{f_{\text{max}}}$$ (6) and $$c_2=\frac{(I_{\text{max}}-I_{\text{final}})}{f_{\text{max}}}.$$ (7)

This definition may seem awkward and overly complex, but it guarantees that certain non-sigmoidal corner cases are eliminated (see Appendix).

The meaning of the parameters in the double-sigmoidal model mirror those of the regular sigmoidal model. The time point t′_mid1 determines (but does not exactly represent) the time at which the intensity has risen to half of its maximum, and the time point t′_mid2 (enforced to be larger than t′_mid1) determines (but does not exactly represent) the time at which the intensity has decayed halfway from its maximum to its final value. The prime in t′_mid1 and t′_mid2 indicates that these time points are not the final parameters of interest. The true midpoints t_mid1 and t_mid2 are calculated numerically by sicegar and reported. Similarly, the parameters a′₁ and a′₂ determine the speed of growth and speed of decay, respectively, but do not represent the slopes at times t_mid1 and t_mid2. The true midpoint slopes are calculated by numerically differentiating f_dsig at each midpoint time. As stated before, the parameters I_max and I_final represent the maximum intensity and the final intensity after decay, respectively.

In addition to the parameters that uniquely determine the sigmoidal and double-sigmoidal functions, the sicegar package calculates other useful metrics that describe aspects of the curves (Fig. 1). For the sigmoidal model, we calculate the starting time of growth by finding the point that a line which passes though the midpoint t_mid with the slope a₁ intersects with the x-axis. This value gives the length of time between the start of observations and the first observable increase in the response variable. Using the same line which passes through t_mid with slope a₁, we can also determine the amount of time that has elapsed during exponential growth, by looking at the time point where the line intersects with the maximum intensity (I_max) and subtracting the starting time. For the double-sigmoidal function the time until decay can be calculated in a similar way (Fig. 1B). We chose to calculate these specific start and end time points for their simple geometric interpretation. They tend to provide a conservative estimate for the total duration of growth or decay.

Model fitting and categorization

We implemented R code that can fit data to the sigmoidal and double-sigmoidal curves defined by Eqs. (1) and (5). Before fitting, data is internally normalized to a [0,1] interval on both the x- and y-axis. The sigmoidal and double-sigmoidal models are then fit to the internally normalized data via likelihood maximization, using the function nls.lm in the minpack.lm package (Elzhov et al., 2015). The minpack.lm package uses a modified Levenberg–Marquardt algorithm (Moré, 1978) to minimize the sum of the squared residuals. To guarantee robust fitting, we run nls.lm multiple times using different starting estimates, chosen at random from a uniform distribution covering the allowed interval of parameters. For each fit, the quality of fit based on Akaike information criterion (AIC) (Akaike, 2011; Burnham & Anderson, 2004) scores is recorded, and the parameter estimates that produce the overall best fit across replicates are reported as the final estimates.

The simplest way to use sicegar is via the function fitAndCategorize, which attempts to both identify the best-fitting model and determine whether this model is sigmoidal, double-sigmoidal, or neither. This function first checks whether there is any signal in the data, by testing whether the maximum observed intensity exceeds a user-defined threshold and the data range also exceeds a user-defined threshold. If either test fails, the data are labeled as “no signal.”

Next, the function fits both a sigmoidal and a double-sigmoidal curve to the data and then assess which provides the better fit. This process considers a number of different factors. First, it assesses whether both models could actually be fitted (i.e., the fitting did not fail). Second, the AIC scores of the fitted models are examined to make sure they are smaller than a threshold value (the default is −10). Third, the starting time points of each model are checked to ensure that they are positive and the starting intensities are examined to ensure that they are smaller than a threshold (the default is 0.05). Fourth, the ratio of the models’ intensity prediction at the last observation to the models’ maximum intensity prediction is tested to ensure that it is over 0.85 in the case of the sigmoidal model and below 0.75 in the case of the double-sigmoidal model. These various checks and cutoffs are heuristics which we developed by manually inspecting thousands of fits in the context of our single-cell virology study (Guo et al., 2017). Some of these (in particular the AIC cutoff) may not be appropriate for a very different dataset and may require adjustment. Therefore, none of these cutoffs are hardcoded into sicegar and they can all be adjusted as needed.

If any of the above tests fail for either the sigmoidal or double-sigmoidal fit, then the data is not categorized as being described by that model. If both models fail any of these tests the data is categorized as being “ambiguous.” If only one model passes these tests while the other fails, then the data is considered to be best described by that model. If both of the fits pass all of these test a decision is made about which model best describes the data based on their AIC scores. If the sigmoidal AIC score is smaller than the double-sigmoidal AIC score then the data is labeled as “sigmoidal,” if it is larger then the data is labeled as “double-sigmoidal.” There is also an option for favoring one or the other model by adding a bonus value to the AIC score of either the sigmoidal or the double-sigmoidal model. By default, the bonus is turned off (set to zero).

Users can modify the complete decision process, as well as the specific parameters to adjust the package to their own needs, as explained in detail in the vignette “Identifying the best-fitting model category” supplied with the R package.

Simulations

To assess sicegar performance, we simulated sigmoidal and double-sigmoidal data with different levels of noise. We first generated sets of parameters for both sigmoidal and double-sigmoidal curves. For the sigmoidal model, the maximum and midpoint parameters were chosen uniformly at random from the intervals [0.3, 20] and [3, 27], respectively. Slopes were generated by first drawing a random slope angle θ uniformly from the interval [0, π/2] and then converting into a slope via tan(θ). For the double-sigmoidal model, the maximum, the first midpoint, the decreasing slope, the distance between midpoints, and the final asymptote intensity ratio parameters were chosen uniformly at random from the intervals [0.3, 20], [3, 26], [0.001, 40], [1, 27—first midpoint], and [0, 0.85], respectively. The increasing slope was generated as in the sigmoidal model, by uniformly sampling the θ interval [0, π/2]. We excluded all randomly generated parameter sets that by definition would be classified as ambiguous by sicegar, for example because the difference between the maximum and minimum intensity fell below our pre-defined cutoff.

For each parameter set, we simulated data for a time interval of 27 units, from t = 3 to t = 30. We used five different temporal sampling approaches, and in each case generated time points and then the associated intensity values according to the chosen model and model parameters. The sampling approaches we considered were: (i) equidistant sampling, where we generated 55 time points in increments of 0.5; (ii) random sampling from a uniform distribution, where we generated 55 time points uniformly distributed from 3 to 30; (iii)–(v) random sampling from beta distributions, creating higher sampling densities at the beginning, center, and end of the time course. Here, we set the α and β parameters of the beta distributions to {α = 0.5, β = 1.5} (mode at the beginning), {α = 2, β = 2} (mode at the center), and {α = 1.5, β = 0.5} (mode at the end), respectively, and sampled 55 time points from the respective beta distribution and then multiplied them by 27 and added 3, so all time points fell into the interval from 3 to 30.

For each generated intensity curve, we then created noisy variants. We used 11 distinct noise amplitudes, ranging from 0% to 150% of the maximum value of the initial simulated data. To generate additive noise, we created vectors of 55 random numbers chosen uniformly from −0.5 to 0.5. We then multiplied each random number with the noise amplitude (number between 0 and 1.5) and the maximum intensity of the dataset. Finally, we added these noise values to the intensity values of the initial samples and obtained noisy samples. To generate multiplicative noise, we generated vectors of 55 random numbers 2^x, where x was uniformly distributed between −1 and 1. We then multiplied each random number with the noise amplitude, and we then multiplied the noise vector element-wise with the original intensity vector.

We generated 50 replicate parameter sets for each combination of model (sigmoidal or double-sigmoidal), temporal sampling approach, noise amplitude, and noise type (additive or multiplicative). We generated a new set of random parameters in each case, so that in total we evaluated 50 × 2 × 5 × 11 × 2 = 11,000 different parameter sets.

For each simulated sample, we used sicegar to estimate both the shape of the curve (“sigmoidal,” “double-sigmoidal,” or “ambiguous”) and the curve’s original parameter values. To estimate the algorithm’s accuracy in fitting a model to the data, we calculated the mean absolute error between the initial (without noise) and predicted intensity vectors, and we normalize this number relative to the maximum intensity of the initial sample. We calculated the normalized mean absolute error regardless of whether the algorithm predicted the category of the sample correctly or not.

Package availability and documentation

The sicegar package is implemented in R and available from CRAN at https://CRAN.R-project.org/package=sicegar. The sicegar source code is available on github at https://github.com/wilkelab/sicegar. Extensive user documentation is available in multiple vignettes distributed as part of the package.

The package depends on existing R packages dplyr (Wickham & Francois, 2015), minpack.lm (Elzhov et al., 2015), fBasics (Wuertz, Setz & Chalabi, 2014), and ggplot2 (Wickham, 2009).

Performance against simulated data

Results from our simulation benchmark are shown in Fig. 2 (for equidistant temporal sampling and additive noise only) and in Fig. S1 (for all five sampling schemes and both additive and multiplicative noise). We first asked how well sicegar recovers the type of the curve (sigmoidal or double-sigmoidal) at varying noise levels (Fig. 2A; Fig. S1A). For equidistant sampling with additive noise, we found that curves that were simulated with sigmoidal input were virtually never classified as double-sigmoidal (Fig. 2A). However, starting at noise levels of around 45% (corresponding to a noise amplitude of half the total amplitude of the sigmoidal curve), sicegar classified a fraction of the curves as ambiguous. The fraction of curves classified as ambiguous increased with increasing noise level until virtually all curves were classified as ambiguous at noise levels of around 100%. Results were similar for curves simulated as double-sigmoidal, except that a small fraction of curves was always identified as sigmoidal, regardless of the noise level (Fig. 2A). These were curves for which the parameter values happened to cause very little decay at large times, either because the final asymptotic value was close to the maximum or because decay started late relative to the maximum time point simulated. Finally, at high noise levels, double-sigmoidal curves can appear sigmoidal, and sicegar classified those curves as such. In general, decreased performance with increased noise level is an expected behavior (Spiess & Neumeyer, 2010), but our algorithm’s tendency to label misclassified cases as ambiguous instead of making wrong decisions shows the robustness of the package’s decision process.

Results were broadly similar for other temporal sampling regimes and multiplicative instead of additive noise (Fig. S1A), with two caveats. First, sicegar was generally more robust to multiplicative noise than to additive noise, though some sigmoidal curves were classified as double-sigmoidal rather than ambiguous under multiplicative noise. Second, for non-equidistant temporal sampling, the performance was best when sampling was uniform or biased toward the middle of the time interval, rather than toward the beginning or end of the interval.

We also compared the fits sicegar obtained from noisy data to the initial curves before applying noise. For each fit, we calculated the normalized mean absolute error, which is the mean of the absolute differences of fitted and original intensity levels, normalized by the maximum original intensity for all five different sampling regimes and for two different noise regimes. We found that these errors were near zero for low noise levels and increased gradually as noise increased (Fig. 2B; Fig. S1B). This was true for both additive and multiplicative noise and for all five temporal sampling regimes. Overall, we can conclude that our algorithm results in reliable fits, that it fails gradually with increasing noise levels, and that it is conservative in assessing whether it has correctly identified a sigmoidal or double-sigmoidal curve or not.

Comparison to existing packages

On the face of it, fitting a sigmoidal curve to time-course data appears to be a simple and solved problem. The simplest sigmoidal curve is determined by only three parameters, and many non-linear curve-fitting algorithms exist that can be used to fit a simple three-parameter function. However, when we attempted to fit thousands of curves to experimentally measured data, we found that achieving numerical reliability was non-trivial. While a simple algorithm will work most of the time, it will fail eventually on cases that by manual inspection have a clear solution. The problem was further exacerbated for double-sigmoidal curves, where even guaranteeing that the curve is double-sigmoidal at all times was non-trivial (see Appendix). Finally, for our original application of single-cell poliovirus infections (Guo et al., 2017), we wanted to classify sigmoidal and double-sigmoidal curves in an automated fashion, and this problem goes beyond just fitting the curves to requiring a heuristic that can automatically classify curves and also determine when the classification is ambiguous.

We wrote the sicegar package to solve these challenges. However, sicegar is not the only package with the ability to fit data to such non-linear models. The qpcR (Ritz & Spiess, 2008) package fits linear, exponential, a number of sigmoidal, and log–logistic functions to qPCR (realtime quantitative polymerase chain reaction) data. Unlike sicegar this package uses a log–logistic function rather than a double-sigmoidal function to describe data that includes both a growth and a decay phase. A similar package, grofit (Kahm et al., 2010), uses a number of parametric growth curves and a model-free spline approach to estimate parameters from dose to response data. This package does not include a parametric model that would account for a decay phase in the data, though an interface to construct custom models is provided. The intended use of this package is to model pharma-/toxicological data and it produces estimates of values such as half maximum effective concentration (EC50). While sicegar can estimate the midpoint, which corresponds to the EC50, sicegar cannot currently estimate time points corresponding to other specific intensity values along the curve. Another package, FlexParamCurve (Oswald et al., 2012), fits in excess of 30 models to data. While this package is useful for teasing apart exactly which variant of a model best describe dynamics, sicegar is better suited for classifying data as either containing sigmoidal or double-sigmoidal dynamics.

We can think of several future extensions to the sicegar package. First, instead of the three-parameter logistic function, we could describe the sigmoidal model via the five-parameter Richards’ equation (Richards, 1959), which is currently available in the growthmodels package. For completeness, however, we would then also have to develop a double-sigmoidal version of this model. Second, future work could add the functionality required to calculate the time points corresponding to arbitrary percentages of the maximum effective concentration, such as the 5th percent (EC5) and the 95th percent (EC95), which could be used as alternative to the start and end point we currently calculate.

Properties of the sigmoidal function

The sicegar package uses for its sigmoidal fit the Fermi function (Eq. 1). Here we describe the properties of a more general variant of this function, $$I(t)=\frac{I_{\text{max}}-I_{\text{init}}}{1+\exp \left[-a_1(t-t_{\text{mid}})\right]}+I_{\text{init}}.$$ (8) The Fermi function is recovered by setting I_init = 0.

The more general function discussed here has four parameters, each of which has a clear meaning in terms of the shape of the function:

1. I_init represents the initial value of the sigmoidal curve at minus infinity, $$I_{\text{init}}=\underset{t\to -\infty }{\lim }I(t).$$ (9)

2. I_max represents the maximum value of the curve at plus infinity, $$I_{\text{max}}=\underset{t\to \infty }{\lim }I(t).$$ (10)

3. t_mid represents the time point at which the slope of I(t) is maximal, $$t_{\text{mid}}=\text{argmax}\frac{d}{dt}I(t).$$ (11)

4. The parameter a₁, constrained to be positive at all times, is related to the maximum slope of I(t), $$\frac{a_1}{4}(I_{\text{max}}-I_{\text{init}})=\frac{d}{dt}I(t){|}_{t=t_{\text{mid}}}$$ (12)

5. In the following, we will also need the inverse sigmoidal function, $$I(t)=\frac{I_{\text{max}}-I_{\text{final}}}{1+\exp \left[a_2(t-t_{\text{mid}})\right]}+I_{\text{final}},$$ (13)

Properties of the double-sigmoidal function

The double-sigmoidal function is defined in Eq. (5) in the main text. It rises from zero to a maximum value and then decays toward an asymptotic final value, and it is defined for all t. The function is continuous with a continuous first derivative and a discontinuous second derivative.

The parameters I_max and I_final represent the maximum and asymptotic values of the curve, respectively. However, the parameters t′_mid1, t′_mid2, a′₁, and a′₂ do not directly represent midpoints and slopes. Instead, the two midpoints and two slopes need to be calculated numerically.

Derivation of the double-sigmoidal function

The double-sigmoidal as defined in Eq. (5) may seem overly complex and awkward. However, we had to define it in this way to avoid problems that arise with other formulations.

In particular, a seemingly more straightforward double-sigmoidal could be obtained by multiplying together a regular sigmoidal (Eq. 8) and an inverse sigmoidal (Eq. 13): $$I(t)=\left(\frac{I_{\text{max}}-I_{\text{init}}}{e^{-a_1(t-t_{\text{mid1}})}+1}+I_{\text{init}}\right)\left(\frac{1-I_{\text{final}}}{e^{a_2(t-t_{\text{mid2}})}+1}+I_{\text{final}}\right)$$ (14) Even though this function behaves appropriately for most parameter choices, there are parameter regions in which this function is not a double-sigmoidal. This occurs in general when the interval between t_mid1 and t_mid2 is very small, and one of the slopes is very steep while the other one is not. As an example, consider the case shown in Fig. 3, where Eq. (14) generates a local minimum in addition to a local maximum. These kinds of corner cases will be discovered by the fitting algorithm when used on sufficiently many diverse datasets, and therefore we ruled out Eq. (14) as an appropriate choice for our double-sigmoidal function.

However, Eq. (14) will never have a local minimum if we set both I_init and I_final to zero (see next section for proof). Thus, we define the base double-sigmoidal as $$f_{\text{dsig-base}}(t)=\frac{1}{\left(e^{-a_1(t-t_{\text{mid1}})}+1\right)\left(e^{a_2(t-t_{\text{mid2}})}+1\right)}$$ (15) This function has the correct shape in principle, except that the intensity always decays back to zero for large t.

One possible solution to the problem that the function approaches zero for $t\to \infty$ is to cut it into two parts, rescale them separately, and then merge them together again. We cut the function at its maximum, where t = t^*. Thus, the part from t = −∞ to t = t^* represents sigmoidal growth and the part from t = t^* to t = ∞ represents sigmoidal decay. Both parts are rescaled such that the desired I_max and I_final are obtained. This procedure leads to the definition of the double-sigmoidal we use, as defined in Eq. (5). This function has the desired limiting properties and is guaranteed to have only one maximum. However, the cutting and rescaling procedure causes its second derivative to be discontinuous at t = t^*.

The base double-sigmoidal has exactly one local maximum

For the base double-sigmoidal as defined in Eq. (15) we can prove the following Lemma.

By definition, t_mid2 > t_mid1, so without loss of generality we can write t_mid2 = t_mid1 + L, where L is a positive number. Then we obtain $$\frac{d}{dt}\left[\frac{1}{(e^{-a_1(t-t_{\text{mid1}})}+1)(e^{a_2(t-L-t_{\text{mid1}})}+1)}\right]{|}_{t=t^{*}}=0$$ (17)

We define new variables u ≡ t − t_mid1 and u^* ≡ t^* − t_mid1. We also define the function g(u) = f_dsig-base(u + t_mid1). Finding the roots of g(u) is equivalent to finding the roots of f_dsig-base(t), and it involves solving the following equation for u^*: $$\frac{d}{du}g(u){|}_{u=u^{*}}=\frac{d}{du}\left[\frac{1}{(e^{-a_{1}u}+1)(e^{a_2(u-L)}+1)}\right]{|}_{u=u^{*}}=0$$ (18)

After taking the derivative, we can rewrite this equation as $$\frac{e^{(a_1{u}^{*}+a_2(L-u^{*}))}\left[a_1(e^{a_2(L-u^{*})}+1)-a_2(e^{a_1{u}^{*}}+1)\right]}{{(e^{a_1{u}^{*}}+1)}^2{(e^{a_2(L-u^{*})}+1)}^2}=0$$ (19)

For the left-hand side to be equal to zero, the term in the square brackets in the numerator must be equal to zero, because all other terms in the expression are always positive. We next define the term in the square brackets as a function of u, $$h(u)=a_1(e^{a_2(L-u)}+1)-a_2(e^{a_{1}u}+1).$$ (20) The number of roots of h(u) is equal to the number of local extrema of f_dsig-base(t).

To determine the number of roots of h(u), we first note the following limiting properties: $$\underset{u\to -\infty }{\lim }h(u)=\infty$$ (21) $$\underset{u\to \infty }{\lim }h(u)=-\infty$$ (22)

According to the mean value theorem, h(u) must thus have at least one root. Next, we note that h(u) is strictly decreasing: $$\frac{d}{du}h(u)=-a_1{a}_2{e}^{a_2(L-u)}-a_1{a}_2{e}^{a_{1}u}<0$$ (23)

Since the derivative of h(u) is not equal to zero for any u we can use Rolle’s Theorem to conclude that the function cannot have more than one root. In conclusion, h(u) has exactly one root, which implies that f_dsig-base(u) has exactly one local extremum.

Finally, we demonstrate that the local extremum of f_dsig-base(t) is a local maximum, by inspecting the second derivative of g(u) at u = u^*.

The first term in the square brackets is equal to h(u)² and thus is zero at u = u^*. The two remaining terms inside the square brackets are strictly negative, and the term outside the square brackets is strictly positive. Therefore, g(u^*) < 0, which proves that the local extremum is a local maximum.