Parámetros de modelos de crecimiento poblacional

Vamos a estudiar dos modelos no-lineales que se usan para describir el crecimiento poblacional de organismos, y vamos a usar el ejemplo de crecimiento bacteriano.

Función Gompertz

La forma general de la función Gompertz, es la siguiente: \[y = ae^{-bc^x}\]

Crecimiento Logístico

Esta es la forma general de la función de crecimiento logístico:

\[N_t = \frac{K}{1 + (\frac{K - N_0}{N_0})e^{-rt}}\]

Ambas son no-lineales, y antes debemos encontrar valores iniciales de los parámetros.

Datos

Utilizaremos los datos de crecimiento de un cultivo de Salmonella sp. (ver proyecto de modelo no-lineal).

library(readxl)
gbac <- read.csv("bacteria19.csv")

Gráfica

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

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

Una forma de 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 más fácilmente, primero creamos una función R para la ecuación Gompertz.

Gompertz <- function(x, y0, ymax, k, lag){
      result <- y0 + (ymax -y0)*exp(-exp(k*(lag-x)/(ymax-y0) + 1) )
      return(result)
}

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

y0 = 4.5
ymax = 10.2
lag = 5
k = 0.2

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

# modelo
Gomp1 <-  nls(Logc ~ Gompertz(Time, y0, ymax, k, lag),
               data = gbac,
               start = list(y0=y0, ymax=ymax, k=k, lag=lag))

# parámetros y estadisticas
summary(Gomp1)
## 
## Formula: Logc ~ Gompertz(Time, y0, ymax, k, lag)
## 
## Parameters:
##      Estimate Std. Error t value Pr(>|t|)    
## y0    4.43112    0.07977   55.55  < 2e-16 ***
## ymax 10.35959    0.06175  167.78  < 2e-16 ***
## k     0.87926    0.03772   23.31  < 2e-16 ***
## lag   7.66630    0.50660   15.13 8.96e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1736 on 24 degrees of freedom
## 
## Number of iterations to convergence: 10 
## Achieved convergence tolerance: 6.197e-06

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, 100, by=0.1)
pred <- predict(Gomp1, list(x=time.levels) )
plot(Logc ~ Time, data=gbac, xlim=c(0,100), ylim=c(0, 12),
     xlab="time", ylab="logCFU")
lines(gbac$Time, pred, 
      lty = 1, col = "red")

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

library(growthcurver)

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

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

model.bac <- SummarizeGrowth(gbac$Time, gbac$Logc, bg_correct = "none")
model.bac$vals
## k    k_se    k_p n0  n0_se   n0_p
## 10.534   0.171   7e-29   3.491   0.179   1e-16
## 
## r    r_se    r_p sigma   df  t_mid
## 0.094    0.007   9e-13   0.437   25  7.497
## 
## t_gen    auc_l   auc_e
## 7.406    922.81  915.104
model.bac$model
## Nonlinear regression model
##   model: n ~ k/(1 + ((k - n0)/n0) * exp(-r * t))
##    data: d
##        k       n0        r 
## 10.53397  3.49142  0.09359 
##  residual sum-of-squares: 4.769
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 1.49e-08

Significado e interpretación de los parámetros

A continuación la lista de parámetros poblacionales que se extraen con el paquete growthcurver y que tienen un significado en el análisis del crecimiento de microorganismos:

  • k es la capacidad de acarreo (\(K\)) de la población en el medio y condiciones de crecimiento. k_se es el error estándar y k_p el nivel de significancia del parámetro (la probabilidad de error al rechazar que el parámetro es igual a 0).
  • n0 es el tamaño de la población inicial de bacterias (\(N_0\)). n0_se es el error estándar y n0_p el nivel de significancia del parámetro (la probabilidad de error al rechazar que el parámetro es igual a 0).
  • r es la tasa de crecimiento poblacional (\(r\)). r_se es el error estándar y r_p el nivel de significancia del parámetro (la probabilidad de error al rechazar que el parámetro es igual a 0). No necesariamente igual a la k del modelo Gompertz.
  • sigma es el error estándar residual, una medida de la bondad de ajuste de los parámetros del modelo a los datos (mejor cuanto más pequeño sea).
  • df son los grados de libertad, calculados como el número de datos (\(n\)) menos el número de parámetros del modelo (\(n - 3\)).
  • t_mid es el tiempo cuando la población alcanza la mitad de la capacidad de acarreo (\(\frac{1}{2}K\)).
  • t_gen es el tiempo de generación más rápida o tiempo de duplicación de la población. Se puede calcular como \(\frac{ln2}{r}\).
  • auc_l es el área bajo la curva del modelo logístico y sirve de medida integradora de \(K, r\) y \(N_0\).
  • auc_e es el área bajo la curva empírica, es decir usando los datos medidos del crecimiento bacteriano.

Gráfica de puntos y del modelo

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