1 Where we left off

In the previous laboratory we arrived at the exponential growth formula:

\[ \frac{dN}{dt} = r N \quad\Longrightarrow\quad N(t) = N_0\, e^{rt} \]

With \(N_0 = 100\) water lilies and \(r = 0.1\)/day, let’s see what the model predicts over the long run:

N0 <- 100
r  <- 0.1

cat("In 100 days:", format(N0 * exp(r * 100),  big.mark = ","), "water lilies\n")
## In 100 days: 2,202,647 water lilies
cat("In 200 days:", format(N0 * exp(r * 200),  big.mark = ","), "water lilies\n")
## In 200 days: 48,516,519,541 water lilies
cat("In 365 days:", format(N0 * exp(r * 365),  big.mark = ","), "water lilies\n")
## In 365 days: 7.108019e+17 water lilies

It is absurd. In less than a year, our small lake is supposed to have more water lilies than there are cells in the human body.

The problem is that the exponential model ignores space, light, nutrients and competition. In this laboratory we will fix this shortcoming.

2 The idea of the carrying capacity \(K\)

Every real environment has a limit to the number of individuals it can sustain. We call this limit the carrying capacity and denote it by \(K\).

For our lake, let’s assume:

K <- 10000   # maximum number of water lilies the lake can hold

A note on the choice: The value \(K = 10{,}000\) is illustrative and was chosen for clarity in the plots. In a real study \(K\) is estimated from data or from the physical features of the habitat (lake surface area, average water lily size, available light, etc.).

2.1 \(K\) is not just a “ceiling”

The naive idea would be: “we let the population grow exponentially up to \(K\), and then we stop it.” But that is not what happens in nature — growth slows down gradually as we approach \(K\), long before we reach it.

Why?

3 A mechanism: why \(K\) reduces the rate

Imagine our lake with the water lilies. When the population is small, each water lily has plenty of space around it and can reproduce freely. But as the population grows, the water lilies in the interior of the cluster are surrounded by other water lilies and struggle to expand. Only the water lilies at the periphery have space. The fuller the lake, the smaller the fraction of water lilies that can actually reproduce.

This is one example of a mechanism. In reality several reasons can act together:

  • Space / light: the water lilies in the interior are shaded.
  • Nutrients: competition for root nutrients intensifies.
  • Pathogens: in dense clusters diseases spread quickly.
  • Toxic by-products: they accumulate as the population becomes denser.

All of these share the same qualitative feature: the per capita reproduction rate decreases as the population approaches \(K\).

So \(K\) is not just a limit we bump into. \(K\) modifies the rate of change itself.

4 From physical intuition to mathematical form

4.1 What we want from the “realized” per capita rate

We define the realized per capita rate of increase as a function \(r_{\text{eff}}(N)\) of the current population. Two natural requirements:

  1. When \(N\) is small, there is plenty of space → \(r_{\text{eff}}(N) \to r\) (the intrinsic rate).
  2. When \(N \to K\), space runs out → \(r_{\text{eff}}(N) \to 0\).

4.2 The simplest function that satisfies these: linear

The simplest function that captures the two extremes is a straight line that falls from \(r\) (at \(N=0\)) to \(0\) (at \(N=K\)):

\[ \boxed{\;r_{\text{eff}}(N) \;=\; r\left(1 - \frac{N}{K}\right)\;} \]

A note on the model choice (important): The linear form is a modeling assumption, not a logical necessity. It is the simplest function that captures the correct qualitative behavior. Other possible forms (we will see them at the end):

  • θ-logistic: \(r_{\text{eff}}(N) = r\bigl(1 - (N/K)^{\theta}\bigr)\)
  • Gompertz: \(r_{\text{eff}}(N) = r\,\ln(K/N)\)

All give similar sigmoidal behavior, but with different geometry.

N_grid <- seq(0, K, length.out = 200)
df_pc <- data.frame(
  N        = rep(N_grid, 2),
  rate     = c(rep(r, length(N_grid)),                # constant r (exponential)
               r * (1 - N_grid/K)),                   # linearly decreasing (logistic)
  Model    = factor(rep(c("Exponential:  r constant",
                          "Logistic: r·(1-N/K)"),
                        each = length(N_grid)))
)

ggplot(df_pc, aes(x = N, y = rate, colour = Model)) +
  geom_line(linewidth = 1.2) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  labs(
    title = "Per capita growth rate as a function of population",
    subtitle = "The essential difference between the exponential and logistic models",
    x = "Population N",
    y = expression(r[eff](N))
  ) +
  theme(legend.position = "bottom")

5 The logistic equation

The total rate of change of the population is:

\[ \frac{dN}{dt} \;=\; r_{\text{eff}}(N) \cdot N \;=\; r\,N\left(1 - \frac{N}{K}\right) \]

This is the logistic equation (Verhulst, 1838).

\[ \boxed{\;\frac{dN}{dt} \;=\; r\,N\left(1 - \frac{N}{K}\right)\;} \]

6 Qualitative analysis before we solve

One of the most useful tools in dynamical systems is to examine the behavior before solving analytically. Very often the solution does not even exist analytically — but the qualitative picture is enough for us.

6.1 Equilibrium points

Where \(\dfrac{dN}{dt} = 0\) the population stays stationary. We solve:

\[ r\,N\left(1 - \frac{N}{K}\right) = 0 \quad\Longrightarrow\quad N^* = 0 \;\text{ or }\; N^* = K \]

  • \(N^* = 0\): empty lake (population extinct).
  • \(N^* = K\): full lake (saturation).

6.2 Stability of the equilibrium points

  • If we start close to \(0\) (e.g. \(N=1\)): \(\frac{dN}{dt} > 0\), the population grows — it moves away from 0. So \(N^*=0\) is unstable.
  • If we start close to \(K\) (e.g. \(N=0.9K\)): \(\frac{dN}{dt} > 0\), the population approaches \(K\).
  • If we start above \(K\) (e.g. \(N=1.5K\)): \(\frac{dN}{dt} < 0\), the population decreases toward \(K\).

So \(N^*=K\) is stable — it is an attractor of the system.

6.3 Phase portrait: \(\dfrac{dN}{dt}\) against \(N\)

N_grid <- seq(-500, K * 1.4, length.out = 300)
dN_dt  <- r * N_grid * (1 - N_grid/K)
df_phase <- data.frame(N = N_grid, dN = dN_dt)

# Equilibrium points
eq_unstable <- data.frame(N = 0, dN = 0)
eq_stable   <- data.frame(N = K, dN = 0)

# Flow arrows on the N axis: the sign of (xend - x) gives the direction of motion
arrows_df <- data.frame(
  x    = c(K/4,        3*K/4,        1.2*K),
  xend = c(K/4 + 800,  3*K/4 + 800,  1.2*K - 800),
  y    = 0, yend = 0
)

ggplot(df_phase, aes(x = N, y = dN)) +
  geom_hline(yintercept = 0, colour = "grey50") +
  geom_line(linewidth = 1, colour = "darkgreen") +
  geom_point(data = eq_unstable, aes(x = N, y = dN),
             shape = 21, fill = "white", colour = "darkgreen",
             size = 4, stroke = 1.2) +
  geom_point(data = eq_stable, aes(x = N, y = dN),
             colour = "darkgreen", size = 4) +
  geom_segment(data = arrows_df,
               aes(x = x, xend = xend, y = y, yend = yend),
               arrow = arrow(length = unit(0.25, "cm")),
               colour = "tomato", linewidth = 0.9,
               inherit.aes = FALSE) +
  annotate("text", x = 0,    y = -150, label = "N* = 0\n(unstable)",  vjust = 1) +
  annotate("text", x = K,    y = -150, label = "N* = K\n(stable)", vjust = 1) +
  annotate("text", x = K/2,  y = max(dN_dt) + 30,
           label = "maximum rate\nat N = K/2", vjust = 0) +
  labs(
    title = "Phase portrait of the logistic equation",
    subtitle = "The arrows show the direction of change of N",
    x = "N", y = "dN/dt = r"
  )

The parabola represents the whole dynamics of the system in a single picture:

  • Where the curve is above the axis → \(\frac{dN}{dt} > 0\) → the population increases.
  • Where the curve is below the axis → the population decreases.
  • At the roots we have equilibrium points.
  • The peak of the parabola (at \(N = K/2\)) is the maximum growth rate — that is where the population grows fastest. This is also the inflection point of the time trajectory.

7 Analytical solution (optional for the mathematically inclined)

For those who want to see the solution: we start with separation of variables:

\[ \frac{dN}{N\,(1 - N/K)} = r\,dt \quad\Longleftrightarrow\quad \frac{K\,dN}{N\,(K - N)} = r\,dt \]

Partial fractions:

\[ \frac{K}{N(K-N)} = \frac{1}{N} + \frac{1}{K-N} \]

So:

\[ \int\!\left(\frac{1}{N} + \frac{1}{K-N}\right) dN = \int r\,dt \;\Longleftrightarrow\; \ln N - \ln(K-N) = rt + C \]

\[ \ln\!\left(\frac{N}{K-N}\right) = rt + C \;\Longrightarrow\; \frac{N}{K-N} = A\,e^{rt}, \quad A = \frac{N_0}{K - N_0} \]

Solving for \(N\) we arrive at the closed-form logistic solution:

\[ \boxed{\;N(t) \;=\; \dfrac{K}{1 + \left(\dfrac{K - N_0}{N_0}\right) e^{-rt}}\;} \]

Sanity checks:

  • At \(t = 0\): \(N(0) = K \big/ \left(1 + \tfrac{K-N_0}{N_0}\right) = K \cdot \tfrac{N_0}{K} = N_0\). ✓
  • At \(t \to \infty\): \(e^{-rt} \to 0\), so \(N \to K\). ✓
  • If \(N_0 \ll K\), for small \(t\) the term \(\tfrac{K-N_0}{N_0} \approx \tfrac{K}{N_0}\) is huge, and the behavior is almost exponential (why?). ✓

8 R: numerical solution with deSolve

Now that we know what to expect, let’s simulate:

logistic_model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dN <- r * N * (1 - N / K)
    list(dN)
  })
}

parameters <- c(r = 0.1, K = 10000)
state      <- c(N = 100)
times      <- seq(0, 200, by = 0.5)

out_log <- as.data.frame(
  ode(y = state, times = times, func = logistic_model, parms = parameters)
)

head(out_log)

8.1 Comparison of the analytical and numerical solutions

N0 <- 100
K  <- 10000
r  <- 0.1

logistic_analytical <- function(t, N0, r, K) {
  K / (1 + ((K - N0) / N0) * exp(-r * t))
}

df_log <- data.frame(
  t          = times,
  Analytical = logistic_analytical(times, N0, r, K),
  Numerical  = out_log$N
)

ggplot(df_log, aes(x = t)) +
  geom_line(aes(y = Analytical),  linewidth = 1.2, colour = "steelblue") +
  geom_point(aes(y = Numerical),
             data = df_log[seq(1, nrow(df_log), by = 25), ],
             colour = "tomato", size = 1.8) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey40") +
  annotate("text", x = max(times) * 0.95, y = K * 1.04,
           label = "K = 10,000", colour = "grey40", hjust = 1) +
  labs(
    title = "Logistic growth: analytical (line) vs numerical (dots)",
    x = "Time (days)",
    y = "Population N(t)"
  )

cat("Maximum error:", format(max(abs(df_log$Analytical - df_log$Numerical)),
                             scientific = TRUE), "\n")
## Maximum error: 6.769672e-04

9 The sigmoidal curve and the three stages

Look carefully at the time trajectory. It has three distinct stages:

stages <- data.frame(
  t = times,
  N = logistic_analytical(times, N0, r, K)
)

# Inflection point: N = K/2
t_inflection <- (1/r) * log((K - N0) / N0)

ggplot(stages, aes(x = t, y = N)) +
  geom_line(linewidth = 1.1, colour = "steelblue") +
  geom_hline(yintercept = K,   linetype = "dashed", colour = "grey50") +
  geom_hline(yintercept = K/2, linetype = "dotted", colour = "grey50") +
  geom_vline(xintercept = t_inflection, linetype = "dotted", colour = "tomato") +
  annotate("point", x = t_inflection, y = K/2, size = 3, colour = "tomato") +
  annotate("text",  x = t_inflection + 5, y = K/2,
           label = "Inflection point (N=K/2)\nmaximum rate", hjust = 0) +
  annotate("text", x =  20, y = 1500, label = "1. Almost exponential",     hjust = 0) +
  annotate("text", x =  60, y = 5500, label = "2. Transitional phase\n(acceleration → deceleration)", hjust = 0) +
  annotate("text", x = 130, y = 9300, label = "3. Saturation toward K",       hjust = 0) +
  labs(
    title = "The sigmoidal (S-shaped) logistic curve",
    x = "Time (days)",
    y = "N(t)"
  )

cat("Inflection point at t =", round(t_inflection, 1), "days (where N = K/2 =", K/2, ")\n")
## Inflection point at t = 46 days (where N = K/2 = 5000 )
  1. Almost exponential phase (small \(N\)): the term \(1 - N/K\) is close to 1, so \(\frac{dN}{dt} \approx rN\).
  2. Inflection point (\(N = K/2\)): here the rate \(\frac{dN}{dt}\) is maximized.
  3. Saturation phase (\(N \to K\)): the term \(1 - N/K \to 0\), so the rate decays exponentially toward zero.

10 Comparison: exponential vs logistic

With the same \(r\), the same \(N_0\), but a different \(K\):

df_compare <- data.frame(
  t           = times,
  Exponential = N0 * exp(r * times),
  Logistic    = logistic_analytical(times, N0, r, K)
)

# Long format for ggplot
df_long <- data.frame(
  t       = rep(times, 2),
  N       = c(df_compare$Exponential, df_compare$Logistic),
  Model   = factor(rep(c("Exponential", "Logistic"), each = length(times)))
)

ggplot(df_long, aes(x = t, y = N, colour = Model)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey50") +
  scale_y_log10() +
  annotate("text", x = max(times) * 0.95, y = K * 1.3,
           label = "K = 10,000", colour = "grey40", hjust = 1) +
  labs(
    title = "Exponential vs logistic growth (logarithmic axis)",
    subtitle = "In the first days they are almost identical — then they diverge dramatically",
    x = "Time (days)",
    y = "N(t) [logarithmic scale]"
  ) +
  theme(legend.position = "bottom")

In the first \(\approx 30\) days the two curves are almost indistinguishable. This gives us a warning: in the early stages of a population phenomenon, we cannot easily tell whether it is growing exponentially or logistically — both models fit equally well. The difference emerges only as we approach \(K\).

11 What happens if we start above \(K\)?

times2 <- seq(0, 100, by = 0.5)

scenarios <- data.frame(
  t   = rep(times2, 3),
  N   = c(logistic_analytical(times2, N0 = 100,   r = 0.1, K = 10000),
          logistic_analytical(times2, N0 = 5000,  r = 0.1, K = 10000),
          logistic_analytical(times2, N0 = 15000, r = 0.1, K = 10000)),
  N0  = factor(rep(c("N0 = 100",  "N0 = 5,000",  "N0 = 15,000 (overshoot)"),
                   each = length(times2)))
)

ggplot(scenarios, aes(x = t, y = N, colour = N0)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey50") +
  annotate("text", x = max(times2) * 0.95, y = K * 1.04,
           label = "K = 10,000", colour = "grey40", hjust = 1) +
  labs(
    title = "Logistic growth: three different initial scenarios",
    subtitle = "All converge to K, regardless of the initial condition",
    x = "Time (days)",
    y = "N(t)"
  ) +
  theme(legend.position = "bottom")

Notice that the case \(N_0 = 15,000 > K\) decreases toward \(K\). This follows automatically from the equation: if \(N > K\), then \(1 - N/K < 0\), so \(\frac{dN}{dt} < 0\).

This models, for example, the case where we introduce an excessive number of water lilies into a lake — the population will fall to the carrying level.

12 Estimating parameters from data

In practice biologists have data and want to estimate \(r\) and \(K\). This is done with non-linear regression (nls() in R).

# We simulate data with noise
set.seed(123)
sample_days <- seq(0, 150, by = 10)
true_N      <- logistic_analytical(sample_days, N0 = 100, r = 0.1, K = 10000)
observed_N  <- true_N * exp(rnorm(length(sample_days), 0, 0.08))

# Fitting the logistic model
fit <- nls(observed_N ~ K / (1 + ((K - N0)/N0) * exp(-r * sample_days)),
           start = list(K = 8000, r = 0.05, N0 = 200))

summary(fit)$coefficients
##        Estimate   Std. Error   t value     Pr(>|t|)
## K  1.023666e+04 241.96331612 42.306681 2.599395e-15
## r  9.074180e-02   0.01263879  7.179627 7.158844e-06
## N0 1.601025e+02  91.19055361  1.755692 1.026589e-01
fit_curve <- data.frame(
  t = times,
  N = predict(fit, newdata = data.frame(sample_days = times))
)

ggplot() +
  geom_point(data = data.frame(t = sample_days, N = observed_N),
             aes(x = t, y = N), size = 2.5, colour = "tomato") +
  geom_line(data = fit_curve, aes(x = t, y = N),
            linewidth = 1, colour = "steelblue") +
  geom_hline(yintercept = coef(fit)["K"],
             linetype = "dashed", colour = "grey40") +
  annotate("text", x = max(times) * 0.95,
           y = coef(fit)["K"] * 1.05,
           label = paste0("Estimate K = ", round(coef(fit)["K"])),
           colour = "grey30", hjust = 1) +
  labs(
    title = "Fitting the logistic model to data",
    subtitle = paste0("Estimates: r = ", round(coef(fit)["r"], 3),
                      ", K = ",  round(coef(fit)["K"]),
                      ", N0 = ", round(coef(fit)["N0"])),
    x = "Time (days)",
    y = "Population N(t)"
  )

Practical significance: The nls() estimated \(\hat{r}\) and \(\hat{K}\) come with a standard error — so we have statistical confidence intervals for the parameters. This is crucial for comparing two populations or two experimental conditions.

13 Beyond the logistic: other functional forms

As we mentioned, the linear form \(r_{\text{eff}}(N) = r(1 - N/K)\) is one choice. Two important alternatives:

13.1 θ-logistic

\[ \frac{dN}{dt} = r N \left(1 - \left(\frac{N}{K}\right)^{\theta}\right) \]

  • \(\theta = 1\) → classical logistic.
  • \(\theta > 1\) → the population stays close to \(K\) for a longer time before slowing down (e.g. \(r\)-strategists).
  • \(\theta < 1\) → the slowdown starts early (e.g. species sensitive to competition).

13.2 Gompertz

\[ \frac{dN}{dt} = r N \ln\!\left(\frac{K}{N}\right) \]

Used mainly for the growth of masses (e.g. cancerous tumors, organism development). It slows down earlier than the logistic.

gompertz <- function(t, N0, r, K) {
  K * exp(log(N0/K) * exp(-r * t))
}
theta_logistic <- function(t, N0, r, K, theta) {
  # Numerical solution
  parms <- c(r = r, K = K, theta = theta)
  func  <- function(t, state, parms) {
    with(as.list(c(state, parms)),
         list(r * N * (1 - (N/K)^theta)))
  }
  out <- ode(y = c(N = N0), times = t, func = func, parms = parms)
  as.data.frame(out)$N
}

t_vec <- seq(0, 200, by = 0.5)
df_alt <- data.frame(
  t        = rep(t_vec, 4),
  N        = c(logistic_analytical(t_vec, 100, 0.1, 10000),
               gompertz(t_vec,             100, 0.1, 10000),
               theta_logistic(t_vec,       100, 0.1, 10000, theta = 0.5),
               theta_logistic(t_vec,       100, 0.1, 10000, theta = 3)),
  Model    = factor(rep(c("Logistic (θ=1)",
                          "Gompertz",
                          "θ-logistic (θ=0.5)",
                          "θ-logistic (θ=3)"),
                        each = length(t_vec)),
                    levels = c("Logistic (θ=1)",
                               "θ-logistic (θ=0.5)",
                               "θ-logistic (θ=3)",
                               "Gompertz"))
)

ggplot(df_alt, aes(x = t, y = N, colour = Model)) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = 10000, linetype = "dashed", colour = "grey50") +
  labs(
    title    = "Different models of saturation",
    subtitle = "All models have the same r=0.1, K=10,000, N0=100",
    x = "Time (days)", y = "N(t)"
  ) +
  theme(legend.position = "bottom")

Pedagogical remark: The fact that all these models give different curves is the most important message. On real data the choice of model does not follow automatically — it depends on the biology of the organism and must be checked statistically (e.g. AIC).

14 Exercises

  1. Sensitivity to \(K\). Plot logistic curves with \(r=0.1\), \(N_0=100\) and \(K \in \{1{,}000,\ 5{,}000,\ 10{,}000,\ 50{,}000\}\) on the same graph. What changes in the shape of the curve?

  2. Time to reach \(90\%\) of \(K\). Solve analytically: for which \(t\) does \(N(t) = 0.9 K\) hold? Verify numerically.

  3. Inflection point. Show that the rate \(\frac{dN}{dt}\) is maximized at \(N = K/2\). (Hint: differentiate \(\frac{dN}{dt}\) with respect to \(N\) and find the root.)

  4. From the data to \(r\). You are given measurements:

    day <- c(0, 10, 20, 40, 60, 90, 120, 180)
N   <- c(120, 280, 620, 2400, 5800, 9100, 9750, 9970)

    Fit the logistic model with nls() and compute the confidence intervals for \(\hat{r}\) and \(\hat{K}\) (confint(fit)).

  5. Conceptual. In the discussion we mentioned various mechanisms that can justify the decrease of the per capita rate (space, light, nutrients, pathogens). For the water lilies in our lake, which mechanism do you consider most important and why? How would you verify it experimentally?

  6. Beyond the mean. All the models here are deterministic — there is no chance involved. In real populations, however, there is stochasticity (births/deaths are random events). Try to simulate stochastic logistic growth: at each time step, add to \(\frac{dN}{dt}\) a random term rnorm(1, 0, σ). How does the behavior change for different \(\sigma\)?

15 Next step

So far we have examined a single population. In reality, however, populations interact with each other:

  • Predator–prey: Lotka–Volterra equations
  • Interspecific competition: extended logistic with two or more equations
  • Symbiosis: mutually beneficial interactions

These are systems of differential equations and will be the subject of the next laboratory.

16 R session information

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Athens
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] deSolve_1.42  ggplot2_4.0.3
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.7.3        cli_3.6.6          knitr_1.51         rlang_1.2.0       
##  [5] xfun_0.57          otel_0.2.0         S7_0.2.2           jsonlite_2.0.0    
##  [9] labeling_0.4.3     glue_1.8.1         htmltools_0.5.9    sass_0.4.10       
## [13] scales_1.4.0       rmarkdown_2.31     grid_4.6.1         evaluate_1.0.5    
## [17] jquerylib_0.1.4    fastmap_1.2.0      yaml_2.3.12        lifecycle_1.0.5   
## [21] compiler_4.6.1     RColorBrewer_1.1-3 farver_2.1.2       digest_0.6.39     
## [25] R6_2.6.1           bslib_0.10.0       tools_4.6.1        withr_3.0.2       
## [29] gtable_0.3.6       cachem_1.1.0