The Linear-Softplus and Double-Softplus Growth Models#
The phenotypic library ships two complementary kinetic fits for time-course colony sizing on solid media:
:class:
~phenotypic.analysis.LinearSoftplus— 4 parameters, no saturation ceiling. Use when colonies are still in the linear- growth regime, or when you want the saturation tail discarded as observation noise. Pruning is ON by default.:class:
~phenotypic.analysis.DoubleSoftplus— 4 or 5 parameters, with a softplus saturation ceiling. Use when the data shows a clear carrying-capacity plateau and you want to recoversmax(and optionally fit β).
Both classes share the same lag-phase math, the same σ-resolution machinery, the same Gaussian-prior-on-s0 infrastructure, and the same robust-loss handling. DoubleSoftplus adds a softplus ceiling on top of the shared core and a per-group binary mode dispatch (fitted_beta vs fixed_beta). This page is a conceptual walkthrough covering:
The model equations and what each parameter means biologically.
The two-way mode dispatch on
DoubleSoftplus— how the fitter picks between a 4-parameter clamped fit (fixed β) and a 5-parameter fit with saturation sharpness as a free parameter.Every outlier-robustness mitigation wired into the fits: saturation pruning (
LinearSoftplusonly), robust loss, inverse-variance weighting, the replicate-SEM quantile floor, the \(n=1\) timepoint pooled-std fallback, and the optional Gaussian prior on the inoculum size \(s_0\).
The growth-curve literature has pulled in this direction for years — robust M-estimators for outliers (Huber, 1964), inverse-variance weighting for heteroskedastic timepoints, and weakly-informative priors for ill-posed initial conditions (Gelman et al., 2008). The softplus-lag representation itself and the saturation-clamp form extend the classical 3-phase bacterial-growth framing of Zwietering et al., 1990 with smoothly-differentiable transitions that play nicely with gradient-based optimizers.
1. The Model Equation#
The fit combines a softplus lag with an optional softplus saturation ceiling. The unclamped (growth-only) form is
and the clamped form folds this through a second softplus to impose the carrying-capacity ceiling:
Why softplus rather than a piecewise-linear lag?#
The softplus function
is the canonical \(C^\infty\) approximation of \(\max(x, 0)\): far below \(x=0\) it decays exponentially to \(0\), far above it grows linearly with slope \(1\), and the transition sharpness is set by \(\alpha\). The lag term in the model is a \(v\)-scaled softplus of \((t - \lambda)\), so the curve transitions smoothly from “no growth” to “linear growth at rate \(v\)” without any discontinuous derivatives. That smoothness matters for gradient-based least-squares solvers: a kink in the model creates a kink in the residual surface, and Levenberg–Marquardt / trust-region solvers do not recover gracefully from kinks.
The saturation clamp plays the same role at the ceiling. Subtracting a softplus of the excess \((s_{\max} - s_{\mathrm{lin}})\) from \(s_{\max}\) flattens the curve smoothly once it approaches the carrying capacity. When \(s_{\mathrm{lin}}(t)\) is well below \(s_{\max}\) the clamp is inactive (softplus(large positive) is essentially the identity); when it is well above, the clamp pins the output to \(s_{\max}\) (softplus(large negative) is nearly zero). The transition
sharpness is set by \(\beta\).
Numerical detail#
The implementation uses :func:numpy.logaddexp, the numerically stable form of \(\ln(1 + e^x)\). Naïve evaluation of \(\ln(1 + e^{\alpha(t-\lambda)})\) overflows when \(\alpha(t-\lambda) \gtrsim 700\); logaddexp(0, x) computes the same quantity as \(\max(0, x) + \ln(1 + e^{-|x|})\) and stays well-behaved across the full range the optimizer visits.
2. Parameter Cheat-Sheet#
Symbol |
Role |
Biological interpretation |
Typical scale |
|---|---|---|---|
\(v\) |
Post-lag linear growth rate |
Maximum specific growth rate reached during the exponential/linear phase |
Units of \(s\) per unit \(t\) (e.g. pixels²/hour) |
\(s_0\) |
Intercept of the post-lag line |
Effective inoculum size, back-extrapolated from the linear phase |
Same units as \(s\) |
\(\lambda\) |
Lag duration |
Time at which the curve enters the linear phase |
Same units as \(t\) |
\(\alpha\) |
Lag transition sharpness |
\(\alpha \to \infty\) gives a hard switch-on; \(\alpha \approx 2\) is very gentle |
Dimensionally \(1/t\); fit bounds are \([2, 50]\) |
\(s_{\max}\) |
Carrying capacity |
Plateau imposed by nutrient limitation, space, or toxic by-product accumulation |
Same units as \(s\) |
\(\beta\) |
Saturation transition sharpness |
Same role as \(\alpha\), but for the ceiling |
Dimensionally \(1/s\); fit bounds are \([2, 50]\) in |
Two things worth highlighting:
:math:`s_0` is an intercept, not the :math:`t=0` observation. The softplus lag is not identically zero for \(t < \lambda\) — it is exponentially small but non-zero. So \(s(0) = s_0 + (v/\alpha)\, \ln(1 + e^{-\alpha\lambda})\), which is numerically very close to \(s_0\) for typical \(\alpha\lambda \gg 1\) but not exactly equal. In the limit of a perfectly sharp switch-on \((\alpha \to \infty)\), the softplus collapses to \(\max(t - \lambda, 0)\) and \(s_0\) is recovered exactly as the pre-lag baseline.
:math:`v` and :math:`lambda` are only separately identifiable with enough sub-saturation points. On a fully-clamped curve, many \((v, \lambda)\) pairs produce the same plateau, and the fit can drift into implausible regions. Section 3 discusses how pruning and priors address this.
3. DoubleSoftplus: Per-Group Mode Dispatch#
LinearSoftplus is single-mode by construction — it always fits the 4-parameter unclamped form. DoubleSoftplus is the saturating variant, and within it the fitter picks one of two variants per group, recorded in the DoubleSoftplus_mode column of the results.
Mode |
Trigger |
Free parameters |
|---|---|---|
|
|
\(v, s_0, \lambda, \alpha\) (clamped, \(\beta\) held constant) |
|
|
\(v, s_0, \lambda, \alpha, \beta\) |
When a curve has no clear plateau, it likely belongs to LinearSoftplus instead — using DoubleSoftplus on non-saturating data still produces a fit (with \(\beta\) pinned at the module default and smax set to the per-group observed max), but the inferred carrying capacity will be sampling-window-dependent rather than biologically meaningful.
Shoulder detection#
Shoulder detection lives in _has_saturation_shoulder and is deliberately conservative:
Dynamic-range gate. If the finite-\(y\) range is \(\leq 5\%\) of \(|\mathrm{median}(y)|\), the group is treated as flat noise and no shoulder is declared. This guards against coincidental ratio tests firing on curves that never really grew.
Peak-slope estimation. The gradient \(dy/dt\) is smoothed with a small box filter whose width scales with the series length (
min(5, max(3, len(y)//4))); the maximum of the smoothed gradient becomes the reference “peak slope” — robust to single-point noise spikes that would otherwise inflate a raw max.Terminal-slope collapse. The minimum of the last few raw gradients (
tail_k = max(2, min(3, len(dy_dt) // 5))) is compared againstshoulder_slope_ratio × peak_slope(default \(0.05 \times\) peak). Using the min (not the mean) of the tail keeps the criterion sensitive to short plateaus where just the final \(2\)–\(3\) points have collapsed to near-zero slope.
If the terminal-slope collapse condition fires, the group has a shoulder; otherwise it does not.
When to pick LinearSoftplus vs DoubleSoftplus#
Picking the right class up front is more honest than relying on runtime dispatch to silently switch models. If a fit group has not plateaued within the observed window and you fit a saturated model, pinning \(s_{\max}\) to the observed maximum silently encodes “this is the carrying capacity” — even though the data only supports “growth has not plateaued yet”. Those are very different claims. The split into two classes makes that choice explicit:
Use
LinearSoftpluswhen growth is still in the linear regime. \(s_{\max}\) and \(\beta\) never appear in the output.Use
DoubleSoftpluswhen there is a real plateau worth recovering. The mode column records whether \(\beta\) was fit or held at the default.
4. Outlier Mitigations#
Field colony data is rarely well-behaved. A typical run has:
A handful of catastrophically wrong measurements — bubble artefacts, segmentation errors, contamination spikes, mis-registered timepoints.
Heteroskedastic noise — early timepoints often have a factor-of-10 smaller absolute \(\sigma\) than late ones, yet contribute the same number of residuals in an unweighted fit.
Replicate coverage that varies timepoint-to-timepoint — some timepoints have many plates, others have only one.
An inoculum whose magnitude is known a priori to a few percent but that the optimizer has no access to — so \(s_0\) routinely drifts into biologically implausible regions when the early-lag fit is under-determined.
Both LinearSoftplus and DoubleSoftplus ship the same set of mitigations, with one exception: saturation pruning lives only on LinearSoftplus (in DoubleSoftplus the saturation IS the model — pruning it would defeat the fit). The mitigations are orthogonal axes of the fit, composable, and the library defaults already sit at the robust end of each dial. Several power-user knobs that used to be public __init__ kwargs (stderr_floor_quantile, v_upper, the pruning
slope-ratio and buffer) are now class constants on the shared base — set them by subclass override if you need to deviate.
4.1 Robust Loss (loss, f_scale)#
:func:scipy.optimize.least_squares exposes a family of robust \(\rho\)-functions that down-weight large residuals (Huber, 1964; the scipy reference lists the exact forms). Both LinearSoftplus and DoubleSoftplus default to loss="huber" with f_scale=1.0.
The optimizer minimises \(\tfrac{1}{2}\,\mathrm{cost}\) where
and \(\rho\) is the chosen loss function applied to the squared, rescaled residuals. With \(\rho(z) = z\) (loss="linear") this is the classical least-squares objective; with the Huber \(\rho\), residuals well below the \(f_\mathrm{scale}\) knee contribute like \(z\) (i.e. squared error), while residuals well above it contribute proportionally to \(\sqrt{z}\) (i.e. absolute error, reducing their influence).
f_scale is the knee of the loss. Its units depend on whether the fit is weighted:
Unweighted fit. \(f_\mathrm{scale}\) is measured in the same units as the residual — so \(f_\mathrm{scale} = 1\) means “any residual larger than \(1\) unit of \(s\) starts getting downweighted.”
Weighted fit. Residuals are pre-divided by \(\sigma\) before entering the loss, so \(f_\mathrm{scale}\) is measured in multiples of \(\sigma\) — \(f_\mathrm{scale} = 3\) treats points more than \(\sim\!3\sigma\) off the curve as outliers.
Choice of loss:
"linear"— classical least-squares, no outlier protection."huber"(default) — quadratic inside the knee, linear outside. The least aggressive robust option; recovers least-squares efficiency on clean data while limiting the influence of true outliers."soft_l1"— smooth approximation to \(\ell_1\) loss; similar to Huber but with a smoother transition at the knee."cauchy"— heavy-tailed; strongly suppresses outliers but can introduce multiple local minima."arctan"— bounded influence function (asymptotes to a constant); the most aggressive outlier suppression available in scipy.
4.2 Inverse-Variance Weighting (stderr_label)#
When replicates are available, the fitter divides each residual by a per-timepoint standard error,
so tight-replicate timepoints pull the fit harder than noisy ones. This addresses heteroskedasticity directly: at late timepoints colony variance is often 5–10× the early-timepoint variance, and without weighting those residuals dominate the loss purely because they are larger in absolute terms, not because they are more informative.
Priority order inside _resolve_y_stderr:
User-supplied ``stderr_label``. If the user passes a column name and that column exists in the aggregated group, those \(\sigma\) values are used verbatim. No pool fallback — the user has opted into their own \(\sigma\) semantics.
Auto-derived replicate SEM (when
stderr_label is None). Duringanalyzethe fitter computesdata.groupby(groupby + [time_label])[on].transform('sem')and stores it in a column named{on}_stderr. That column is then read as \(\sigma\) for each aggregated row.No weights (return
None) if neither path is available.
The auto-SEM path is the common case: the user does not need to pre-compute \(\sigma\), and the per-timepoint standard error naturally captures both the measurement noise and any per-plate variability contributing to that timepoint.
4.3 Replicate-SEM Quantile Floor (class constant _STDERR_FLOOR_QUANTILE)#
Pure inverse-variance weighting has a well-known failure mode. When replicates happen to agree by chance — or when the measurement process is discretised so a handful of reps land on identical values — the SEM collapses toward zero and the corresponding \(1/\sigma^2\) weight blows up. A few “lucky” timepoints can then pin the entire fit, because the optimizer will sacrifice a huge amount of fit quality elsewhere to gain a tiny improvement on the \(1/\sigma^2 \to \infty\) residuals.
The shared base class sets _STDERR_FLOOR_QUANTILE = 0.25 (a class constant, not an __init__ kwarg). The procedure is:
Collect all finite, positive \(\sigma\) values in the group.
Compute the \(q\)-th quantile (default \(q = 0.25\)) of that set.
Replace any \(\sigma_i\) below the quantile with the quantile value: \(\sigma_i \leftarrow \max(\sigma_i, \sigma_q)\).
The ratio between the stiffest and softest point in a group is therefore capped at \(\sigma_{\max} / \sigma_q\), which at \(q = 0.25\) is at most a factor of about \(16\) under typical empirical \(\sigma\) distributions. That is large enough to preserve legitimate heteroskedasticity (orders-of-magnitude differences across time) but small enough that a coincidentally-agreeing replicate cannot dominate the fit.
Two notes:
Disabling. Subclass
LinearSoftplus(orDoubleSoftplus) and set_STDERR_FLOOR_QUANTILE = Noneto recover raw inverse- variance weighting. This is appropriate when thestderr_labelcolumn is trusted across its full range, or when you are deliberately probing the effect of the floor.All-zero / all-NaN fallback. When a fit group has no finite positive \(\sigma\) at all (every timepoint has \(n=1\), or noise-free synthetic replicates that all coincide), the model bypasses the weighted pathway entirely and falls back to unweighted residuals. The floor is inapplicable in that case by construction.
The floor applies uniformly to both user-supplied stderr_label columns and the auto-derived replicate-SEM column, so the same calibration works across both paths.
4.4 Pooled Point-Level Std for \(n=1\) Timepoints#
A related pathology: if a fit group has a mix of multi-replicate and singleton (\(n=1\)) timepoints, the singletons get \(\sigma = \mathrm{NaN}\) from the SEM calculation. A naïve fix — filling the NaN with a tiny \(\varepsilon\) — would let the singleton point dominate the fit with \(1/\varepsilon^2\) weight, the same pathology the quantile floor addresses but arriving via a different route.
The fitter’s auto-\(\sigma\) path pre-computes a pooled point-level std alongside the SEM column during analyze. Concretely:
{on}_stderr— per-timepoint standard error, viagroupby(groupby + [time_label])[on].transform('sem').{on}_std_pool— per-timepoint std broadcast and then median-reduced across the fit group:std_pt = groupby(groupby + [time_label])[on].transform('std') pool = std_pt.groupby(groupby).transform('median')
So each fit group gets one scalar: the median of its multi-replicate timepoints’ per-timepoint stds.
Inside _resolve_y_stderr, NaN / zero entries in the SEM vector are filled with the pool value rather than \(\varepsilon\). The \(n=1\) rows end up with \(\sigma \approx\) the typical per-timepoint noise level, so their weight is commensurate with — rather than dominant over — the multi-replicate rows.
Accepted trade-off: median-of-transform ensures that timepoints with more replicates contribute proportionally more rows (and hence more weight) to the pool. This tends to enlarge the pool-\(\sigma\) relative to a simple unweighted median of per-timepoint stds, which is the safer direction — pool-filled \(n=1\) rows will never dominate over their better-supported \(n \geq 2\) neighbours.
Three edge cases:
User-supplied ``stderr_label`` path. The pool is not used here; NaN / zero entries fall back to the \(\varepsilon\)-scaled fill. The user has committed to their own \(\sigma\) semantics, so the library doesn’t silently cross-contaminate them with replicate-derived pool values.
Fully-singleton fit group. If every timepoint has \(n = 1\), the pool is itself NaN (std needs \(n \geq 2\)).
_resolve_y_stderrreturnsNoneand the fit runs unweighted — the same fallback as the “no finite positive \(\sigma\)” case.Quantile floor interaction. Pool-filled entries are not included in the set used to compute the floor quantile — the floor is derived from finite, positive, non-filled \(\sigma\) values only — but they are subsequently lifted to the floor via the final
np.maximum. So pool-filled rows inherit both the pool baseline and any floor that is larger than the pool.
4.5 Saturation Pruning (LinearSoftplus.prune_saturated)#
Pruning lives only on LinearSoftplus — DoubleSoftplus models the saturation explicitly and trimming the plateau would defeat the fit.
Post-saturation points cause two problems for a non-saturating model. First, they are often where contamination overruns a colony or where segmentation confuses neighbour-colony merges, pulling the late-time values upward or downward in biologically meaningless ways. Second, even perfectly clean plateau points are uninformative about the lag-phase parameters \((\lambda, \alpha, s_0)\) — yet they contribute an equal share of residuals unless something trims them off.
_prepare_group implements a hybrid heuristic. A point is flagged as saturated when both of the following are true:
Amplitude gate. \(y \;\geq\; \min(y) + 0.90\,\bigl(\max(y) - \min(y)\bigr)\). The point has climbed into the upper \(10\%\) of the observed dynamic range.
Sustained derivative run. Three consecutive smoothed-gradient values below
_PRUNE_SLOPE_RATIO × peak_slope(default \(0.05 \times\) peak). A point is part of a sustained low-slope run, not just an isolated dip.
Both gates must agree before any trimming happens. Only rows past the latest index satisfying both — plus a small _PRUNE_BUFFER class constant (default \(2\)) of rows kept after the cut so the fit still sees some plateau evidence — are removed.
Why both gates?
The amplitude gate is structurally immune to lag-phase noise. During the lag phase \(y \approx s_0\), far below the \(0.9 \times\) amplitude target, so pruning can never fire on lag-phase points no matter how flat the derivative is.
The sustained-run gate rejects transient mid-growth dips. A single low-slope point in the middle of the linear phase (e.g. a measurement where two consecutive timepoints happened to land at similar values) does not satisfy the “three in a row” requirement.
Two additional guards:
Tail-growth short-circuit. If the mean of the tail-window smoothed slope is still above the threshold, pruning returns the group unchanged. Curves that never plateau within the observation window are not trimmed.
Degenerate-input guard. If the smoothed peak slope is non-finite or non-positive (all-NaN \(y\), or a constant series), pruning returns the group unchanged.
Set prune_saturated=False to disable pruning entirely. Subclass LinearSoftplus and override _PRUNE_SLOPE_RATIO / _PRUNE_BUFFER for tuning; these are class constants, not __init__ kwargs.
4.6 Gaussian Prior on \(s_0\) (s0_prior, s0_prior_cv, s0_prior_sigma)#
When the earliest observed timepoint is already well past the lag (e.g. the first capture is at \(t = 12\) hr and the lag was \(\lambda = 4\) hr), the fit is under-determined: \(s_0\) and \(\lambda\) trade off freely along an approximately linear manifold and the optimizer has no principled basis for choosing one point on the manifold over another. A weakly-informative Gaussian prior on \(s_0\) breaks the symmetry without driving the fit toward implausible regions — the standard weakly-informative-prior argument from Gelman et al., 2008.
Implementation as a virtual residual#
The prior is added to the least-squares residual vector as a single extra row:
Because least-squares with Gaussian noise is equivalent to maximum- likelihood under a Gaussian likelihood, appending a Gaussian-shaped virtual residual is equivalent to maximum a posteriori (MAP) estimation with a Gaussian prior on \(s_0\). The optimizer sees no distinction between the prior residual and the data residuals; the prior simply adds \((s_0 - \mu)^2 / \sigma^2\) to the objective, which is exactly the negative log-density of a Gaussian prior up to additive constants.
s0_prior dispatch (by type)#
The s0_prior argument is polymorphic — the type dispatches which \(\mu\)-estimation strategy is used.
Value |
Behaviour |
|---|---|
|
No prior (default). Residual vector is unchanged. |
|
Data-ground on |
|
Data-ground on the named column. \(\mu\) is the median of |
positive |
Scalar \(\mu\), applied uniformly to every fit group. |
anything else |
Raises :class: |
non-positive or non-finite scalar |
Raises :class: |
The bool branch comes before the numeric branch in the underlying helper because bool is a subclass of int in Python — so isinstance(True, int) is True and an unguarded numeric branch would consume the boolean inputs.
Choosing \(\sigma\): s0_prior_cv vs s0_prior_sigma#
The prior \(\sigma\) is set by exactly one of two mutually exclusive knobs:
Argument |
Behaviour |
|---|---|
|
Coefficient of variation: \(\sigma = \mathrm{cv} \times \mu\) |
|
Absolute σ: \(\sigma = \mathrm{sigma}\), independent of \(\mu\) |
If both are set, construction raises :class:ValueError. If neither is set and the prior is engaged, the helper applies CV \(= 0.05\) as a moderately informative default — the same effective default as the legacy s0_prior_factor=0.05. Pick CV when the prior strength should scale with magnitude (typical for biologically positive quantities); pick absolute σ when you have a hard tolerance budget.
Empirical-Bayes pooling (s0_prior_groupby)#
For column-backed priors (s0_prior = True or a string), s0_prior_groupby lets you pool \(\mu\) estimation over a coarser grouping than groupby. For example, if the fit grouping is ["Media", "Strain", "Plate"] but inoculation spread varies per-media and replicates within a media/strain combination are exchangeable, setting s0_prior_groupby=["Media"] computes the prior \(\mu\) as the median of the earliest-timepoint observations within each media, then broadcasts
that one \(\mu\) to every strain/plate within the media.
This is the classical empirical-Bayes move — Efron & Morris, 1975 and the shrinkage literature that followed — for situations where per-fit-group \(\mu\) estimates would be noisy but the group-level average is stable.
Constraints, verified at analyze time:
s0_prior_groupbymust be a strict subset ofgroupby(coarser, not different). Violating this raises :class:ValueError.Columns referenced in
s0_prior_groupbymust exist in the data.s0_prior_groupbymust be non-empty if supplied. PassNone(the default) to fall back to using the full fitgroupbyfor prior-\(\mu\) estimation.
When the prior does (and does not) engage#
The helper’s is_configured property summarises the engagement logic: the prior attempts to apply when either a column label is set (s0_prior=True or s0_prior=<str>) or a scalar direct_mean is set (s0_prior was a positive number). If neither applies, the prior is disabled.
Even when configured, the prior can still skip a group when:
The broadcast \(\mu\) column is missing for the group (the earliest timepoint had no usable column values — e.g. all NaN).
The resulting \(\mu\) or \(\sigma\) is non-finite or non-positive.
In either case, _inoc_stats(group) returns None and the fit runs without the prior residual for that group. Other groups in the same analyze call are unaffected.
5. Putting It Together#
The mitigations are orthogonal axes of the fit. The library defaults already sit at the robust end of every dial:
Argument |
Lives on |
Default |
Addresses |
|---|---|---|---|
|
both |
|
Catastrophic outliers (spikes, mis-segmentation) |
|
both |
|
Knee of the robust loss |
|
|
|
Uninformative plateau points, tail contamination |
|
both |
|
Auto-derives replicate SEM during aggregation |
|
both |
|
Opt-in; off by default because it changes the objective |
|
both |
|
Relative prior strength |
|
both |
|
Absolute prior σ (mutually exclusive with cv) |
|
both |
|
Empirical-Bayes pooling (opt-in) |
|
|
|
Carrying capacity |
|
|
|
Saturation transition sharpness |
|
|
|
Shoulder-detection threshold for fitted_beta vs fixed_beta |
Class constants (subclass override; not __init__ kwargs):
Constant |
Default |
Lives on |
|---|---|---|
|
|
shared base |
|
|
shared base |
|
|
|
|
|
|
The per-group mode dispatch on DoubleSoftplus (fitted_beta / fixed_beta) runs on top of these and needs no configuration in the common case — \(\beta\) is fit when a shoulder is present, held at the default when no shoulder is detected.
References#
Efron, B., & Morris, C. (1975). Data analysis using Stein’s estimator and its generalizations. Journal of the American Statistical Association, 70(350), 311–319. https://doi.org/10.1080/01621459.1975.10479864
Gelman, A., Jakulin, A., Pittau, M. G., & Su, Y.-S. (2008). A weakly informative default prior distribution for logistic and other regression models. The Annals of Applied Statistics, 2(4), 1360–1383. https://doi.org/10.1214/08-AOAS191
Huber, P. J. (1964). Robust estimation of a location parameter. The Annals of Mathematical Statistics, 35(1), 73–101. https://doi.org/10.1214/aoms/1177703732
Zwietering, M. H., Jongenburger, I., Rombouts, F. M., & van ‘t Riet, K. (1990). Modeling of the bacterial growth curve. Applied and Environmental Microbiology, 56(6), 1875–1881. https://doi.org/10.1128/aem.56.6.1875-1881.1990