The model and its parameterisation

This vignette is the mathematical reference for freqTLS: the four-parameter logistic (4PL) thermal death-time curve, the direct CTmax/z parameterisation, the disjoint-bounds asymptote transform, the relative-versus-absolute survival threshold, and the algebraic bridge to bayesTLS. The thermal-load-sensitivity framework was introduced by Daniel W. A. Noble, Pieter A. Arnold, and Patrice Pottier in bayesTLS. The package’s implementation contract is tested against these equations; this vignette mirrors it and checks the bridge identities numerically (no Stan required). If you already understand the matched-configuration bridge to bayesTLS, you can skip straight to the worked examples in vignette("freqTLS") and vignette("comparing-to-bayesTLS").

library(freqTLS)
#> freqTLS 0.1.0
#> Please cite: Noble DWA, Arnold PA, Nakagawa S & Pottier P (2026) A flexible
#>   modelling framework for estimating thermal tolerance and sensitivity.
#>   bioRxiv. doi:10.64898/2026.07.16.738378
#> Run  citation("freqTLS")  for all entries.
#> 
#> Tutorial & online vignette: https://itchyshin.github.io/freqTLS/

The 4PL thermal death-time curve

Let d be the exposure duration and \log d = \log_{10} d. Survival probability follows the descending four-parameter logistic

p = \mathrm{low} + \frac{\mathrm{up} - \mathrm{low}}{1 + \exp\!\big(k\,(\log_{10} d - \mathrm{mid})\big)}, \qquad k = \exp(\ln k).

\mathrm{low} is the lower asymptote (survival at long exposures),
\mathrm{up} is the upper asymptote (survival at short exposures),
k > 0 is the steepness; its internal coordinate is the natural logarithm \ln k, and
\mathrm{mid} is the midpoint on the \log_{10} d axis.

Survival is high at short durations and falls toward \mathrm{low} at long durations.

The direct `CTmax`/`z` parameterisation

freqTLS lets temperature enter the midpoint directly through the two headline quantities. With z_i = \exp(\eta_{\log z, i}) and \mathrm{CTmax}_i a function of the design,

\mathrm{mid}_i = \log_{10}(t_\mathrm{ref}) - \frac{T_i - \mathrm{CTmax}_i}{z_i}.

This parameterisation makes the two biological quantities direct model parameters, so profile likelihood can constrain them directly:

At T_i = \mathrm{CTmax}_i, the midpoint is exactly \log_{10}(t_\mathrm{ref}) — i.e. CTmax is the temperature at which the threshold crossing occurs at the reference time t_\mathrm{ref}.
The slope of the midpoint in temperature is \partial\,\mathrm{mid}/\partial T = -1/z, so z is the change in temperature per order of magnitude (a one-unit change in \log_{10} duration) — the thermal sensitivity in degrees Celsius per 10-fold change in exposure duration.

Because CTmax and z are model parameters (not derived afterwards), each can be profiled directly. We can check the two properties on a fitted model:

set.seed(1)
dat <- simulate_tls(family = "binomial", CTmax = 36, z = 4, seed = 1)
fit <- fit_tls(dat, y = survived, n = total, time = duration, temp = temp,
               family = "binomial", tref = 60)
ct <- get_ctmax(fit)$estimate
zz <- get_z(fit)$estimate

# (1) midpoint at temp = CTmax equals log10(tref) = log10(60)
predict(fit, data.frame(temp = ct, duration = 1), type = "midpoint")
#> [1] 1.778151

# (2) midpoint slope in temperature is -1/z
mids <- predict(fit, data.frame(temp = c(ct, ct + 1), duration = 1), type = "midpoint")
c(slope = diff(mids), minus_one_over_z = -1 / zz)
#>            slope minus_one_over_z 
#>       -0.2501231       -0.2501231

Disjoint-bounds asymptotes

The asymptotes use the bayesTLS compute_4pl_bounds() parameterisation: the feasible band [\ell, u] (default [0, 1]) is split at its midpoint, and low and up each map an unconstrained coefficient onto one half, so \mathrm{up} > \mathrm{low} holds automatically:

\mathrm{low} = \ell_{\min} + w_\mathrm{low}\,\mathrm{logit}^{-1}(\beta_\mathrm{low}), \qquad \mathrm{up} = u_{\min} + w_\mathrm{up}\,\mathrm{logit}^{-1}(\beta_\mathrm{up}).

low is confined to the lower half-band [\ell_{\min}, \ell_{\max}] and up to the upper half-band [u_{\min}, u_{\max}] (the bands meet, with a tiny separating gap, at the midpoint). Any (\beta_\mathrm{low}, \beta_\mathrm{up}) therefore gives a valid ordered pair \ell < \mathrm{low} < \mathrm{up} < u. This is unconstrained and smooth, which keeps the optimiser and the profiles well-behaved, and it matches bayesTLS exactly so the two packages share the asymptote contract.

Under this parameterisation up is fitted through its own coordinate \beta_\mathrm{up}, just as low is fitted through \beta_\mathrm{low}. The limitation concerns interval computation, not fitting: freqTLS does not yet profile \beta_\mathrm{up}, so it reports a delta-method Wald interval for up (or a bootstrap interval when requested). The full table of parameters and links:

Natural parameter	Internal coordinate	Link
`low`	`beta_low`	logit (onto the lower half-band)
`up`	`beta_up`	logit (onto the upper half-band)
`k`	`beta_logk`	natural log (`ln`)
`CTmax` (per group)	`beta_CT[g]`	identity
`z` (per group)	`beta_logz[g]`	log
`phi` (beta-binomial)	`log_phi`	log

Relative versus absolute thresholds

A “lethal time” such as the LT_{50} is the duration at which survival crosses a target probability. There are two conventions. A relative threshold is a position between the fitted asymptotes; an absolute threshold is a fixed survival probability on the response scale:

Relative (the freqTLS default): the target is interpreted relative to the fitted asymptotes, so relative fraction 0.5 means halfway between low and up — which is exactly the 4PL midpoint \mathrm{mid}. This is the configuration that matches the bayesTLS target_surv = "relative" setting and is what the benchmark uses (see vignette("comparing-to-bayesTLS")).
Absolute: the target is an absolute survival probability, requiring p to lie strictly between low and up.

derive_lt() solves the 4PL for the duration at which survival reaches its numeric absolute probability p, which must lie strictly between low and up. At the relative midpoint, use p = (low + up) / 2; plot_tdt_curve() uses that midpoint by default and can plot another valid absolute p. On the \log_{10} d axis the relative-threshold crossing is the line \log_{10} d = \mathrm{mid}(T) = \log_{10}(t_\mathrm{ref}) - (T - \mathrm{CTmax})/z, whose slope is -1/z — the classic log-linear thermal death-time line.

# Absolute 50% survival at three temperatures, in hours
derive_lt(fit, p = 0.5, temp = c(35, 36, 37))
#> [1] 1.6998017 0.9555978 0.5372199

The bridge to bayesTLS

For the matched constant-shape midpoint configuration (temp_effects = "mid"), bayesTLS parameterises the midpoint as a line in temperature,

\mathrm{mid}(T) = b_{\mathrm{mid,Intercept}} + b_{\mathrm{mid},T_c}\,(T - \bar T),

and reads the thermal sensitivity and critical thermal maximum off that line:

z = -\frac{1}{b_{\mathrm{mid},T_c}}, \qquad \mathrm{CTmax}(t_\mathrm{ref}) = \bar T + \frac{\log_{10}(t_\mathrm{ref}) - b_{\mathrm{mid,Intercept}}}{b_{\mathrm{mid},T_c}}.

Expanding the freqTLS midpoint as a line in T gives the inverse map:

\beta_1 = -\frac{1}{z}, \qquad \beta_0 = \log_{10}(t_\mathrm{ref}) + \frac{\mathrm{CTmax} - \bar T}{z},

z = -\frac{1}{\beta_1}, \qquad \mathrm{CTmax} = \bar T + \frac{\log_{10}(t_\mathrm{ref}) - \beta_0}{\beta_1}.

These are the midpoint-line identities used by bayesTLS. The reparameterisation shares (low, up, k) and is smooth and invertible while z is finite. Therefore a matched constant-shape, relative-threshold likelihood has the same fitted curve under either coordinate system. Current bayesTLS also offers direct CTmax/z parameterisation and an absolute-threshold option; numerical comparisons must still match data, formulas, bounds, threshold, and random-effect structure. freqTLS’s distinct contribution is maximum-likelihood fitting with direct profile-likelihood targets.

We can verify the bridge numerically without bayesTLS. Fit the midpoint line directly from the model’s own predicted midpoints (linear in T), then recover CTmax and z from the slope and intercept:

Tbar <- mean(c(35, 37))               # any centring temperature works
mids <- predict(fit, data.frame(temp = c(35, 37), duration = 1), type = "midpoint")
beta1 <- diff(mids) / 2               # slope in T  (= -1/z)
beta0 <- mids[1] - beta1 * (35 - Tbar) # intercept at T = Tbar

z_recovered     <- -1 / beta1
ctmax_recovered <- Tbar + (log10(1) - beta0) / beta1

rbind(
  fitted    = c(CTmax = ct, z = zz),
  recovered = c(CTmax = ctmax_recovered, z = z_recovered)
)
#>              CTmax        z
#> fitted    28.81676 3.998031
#> recovered 35.92586 3.998031