Parámetros de modelos de crecimiento poblacional

Función Gompertz

\[y = ae^{-bc^x}\]

Crecimiento Logístico

\[N_t = \frac{K}{1 + (\frac{K - N_0}{N_0})e^{-rt}}\] Ambas son no-lineales, y algo difíciles para encontrar unos valores iniciales de los parámetros.

Datos

library(readxl)
data <- read_excel("data/ComBaseExport20.xlsx", 
    sheet = "Logcs")
Time <- as.numeric(data$Time)
Logc <- as.numeric(data$Logc)
gbac <- data.frame(Time, Logc)

Gráfica

library(ggplot2)
ggplot(gbac, aes(x = Time, y = Logc)) + geom_point(alpha=0.7)

Método 1: Parámetros del modelo Gompertz modificado

La ecuación re-paremetrizada del modelo Gompertz es la siguiente: \[y = y_0 + (y_{max} - y_0)*exp(-exp(k*(lag - x)/(y_{max} - y_0) + 1))\] Para utilizar el comando nls meas fácilmente, primero creamos una función R para la ecuación Gompertz.

Luego obtenemos unos valores iniciales para los parámetros, a partir de la gráfico y/o los datos.

Finalmente se calculan los parámetros del modelo ajustado usando nls.

# definir función
Gompertz <- function(x, y0, ymax, k, lag){
      result <- y0 + (ymax -y0)*exp(-exp(k*(lag-x)/(ymax-y0) + 1) )
      return(result)
}
# modelo
Gomp1 <-  nls(Logc ~ Gompertz(Time, y0, ymax, k, lag),
               data = gbac,
               start = list(y0=3, ymax=9, k=0.1, lag=4))
# parametros
coefs <- coef(Gomp1)
y0 <- coefs[1]
ymax <- coefs[2]
k <- coefs[3]
lag <- coefs[4]
y0
##      y0 
## 3.06206
ymax
##     ymax 
## 10.95147
k
##         k 
## 0.5136992
lag
##      lag 
## 10.95963

Gráfica de puntos y del modelo

Se grafican los puntos originales y valores de la predicción con la función y los parámetros calculados.

time.levels <- seq(0, 50, by=0.1)
pred <- predict(Gomp1, list(x=time.levels) )
plot(Logc ~ Time, data=gbac, xlim=c(0,50), ylim=c(0, 10),
     xlab="time", ylab="logCFU")
lines(time.levels, Gompertz(time.levels,
                            y0 = y0, ymax = ymax, k = k, lag=lag), 
      lty = 1, col = "red")

Método 2: Paquete especializado para curvas de crecimiento bacteriano con modelo logístico

library(growthcurver)

Función SummarizeGrowth para obtener los parámetros del modelo

Con este método no es necesario estimar valores iniciales, y se puede aplicar a muchas curvas de crecimiento a la vez.

Cálculo de los parámetros del modelo y de análisis de crecimiento bacteriano.

model.bac <- SummarizeGrowth(gbac$Time, gbac$Logc, bg_correct = "none")
model.bac$vals
## k    k_se    k_p n0  n0_se   n0_p
## 13.539   6.259   8e-02   2.612   0.426   2e-03
## 
## r    r_se    r_p sigma   df  t_mid
## 0.046    0.02    7e-02   0.706   5   31.439
## 
## t_gen    auc_l   auc_e
## 15.229   293.846 294.565
predict(model.bac$model)
## [1] 2.612287 3.017663 3.465991 3.955883 4.484192 5.784473 7.922879 9.470132
model.bac
## Fit data to K / (1 + ((K - N0) / N0) * exp(-r * t)): 
##      K   N0  r
##   val:   13.539  2.612   0.046
##   Residual standard error: 0.7060364 on 5 degrees of freedom
## 
## Other useful metrics:
##   DT 1 / DT  auc_l   auc_e
##   15.23  6.6e-02 293.85  294.56

Gráfica de puntos y del modelo

plot(model.bac, xlab="time", ylab="logCFU")