Modelos Empíricos y Mecanicistas

Los modelos estadísticos (empíricos) suelen estar basados en regresiones. Ellos proporcionar un resumen cuantitativo de las relaciones observadas entre un conjunto de variables medidas.

Un modelo mecanicista (determinístico) comienza con una descripción de cómo la naturaleza podría funcionar, y procede de esta descripción a un conjunto de predicciones que relacionan las variables independientes y dependientes.

Pasos y Herramientas para Formular Modelos Mecanicistas en Sistemas Biológicos y realizar Simulaciones

Sistemas que Cambian con el Tiempo (o el Espacio)

El estudio de procesos biológicos es el estudio de cambios que ocurren en el tiempo y en el espacio, por lo tanto los modelos que usamos son modelos dinámicos. Por esto, en esta parte del curso estaremos trabajando principalmente con ecuaciones diferenciales, en modelos continuos. Más adelante haremos algunos ejercicios con ecuaciones de diferencia y matrices, en modelos discretos. Recomiendo revisar el concepto en Khan Academy.

Modelo de Crecimiento Poblacional

Digamos que queremos estudiar el crecimiento de una población de bacterias u otro organismos de reproducción continua (asexual preferiblemente). Partimos de una población inicial \(N_0\) y por nuestro conocimientos de cálculo, podemos representar la velocidad instantánea de crecimiento como: \[\frac{dN}{dt}\] Vamos a considerar que la velocidad de crecimiento es una constante \(a\): \[\frac{dN}{dt} = a\] Podemos resolver esta ecuación para expresar el tamaño de la población \(N\) en función del tiempo \(t\) (pueden utilizar Wolfram Alpha :
\[N(t) = at + C_1\] Nuestro valor inicial \(N_0\) es la constante \(C_1\):
\[N(t) = at + N_0\] Vamos a utilizar el paquete deSolve para integrar la ecuación diferencial y observar el comportamiento al graficarla.

library(deSolve)
library(ggplot2)
# tiempo
time <- seq(from=0, to=30, by = 0.1)
# valor inicial de N
N <- 100
# valor del parámetro a, velocidad de crecimiento
par_a <- 2
# función a integrar
crec1 <- function(t, N, par_a){
        with(
            as.list(c(N, par_a)),{
                dx <- par_a
                return(list(dx))
            }
            )
}
# integración
out <- ode(y = N, times = time, func = crec1, parms = par_a)
# cambiando el formato del archivo de salida
ninteg <- as.data.frame(out)
colnames(ninteg) <- c("time","N")
# valores de N
head(ninteg)
##   time     N
## 1  0.0 100.0
## 2  0.1 100.2
## 3  0.2 100.4
## 4  0.3 100.6
## 5  0.4 100.8
## 6  0.5 101.0
# gráfica
ggplot(data = ninteg, aes(x=time, y=N)) + geom_point()

Comparando este resultado con la experiencia sobre observaciones de como procede el crecimiento de organismos como las bacterias este modelo no parece representarlo. Vamos ahora a considerar que el crecimiento poblacional depende del valor de \(N\), ya que mientras más individuos hay en la población, ocurrirá más reproducción (feedback); también vamos a considerar que el parámetro \(a\) representa el crecimiento per capita, \(r\): \[\frac{dN}{dt} = rN\] Resolviendo analíticamente la ecuación diferencial, tenemos: \[N(t) = N_0e^{rt}\] Usamos nuevamente el paquete deSolve:

# tiempo
time <- seq(from=0, to=30, by = 0.1)
# valor inicial de N
N <- 100
# valor del parámetro a, velocidad de crecimiento
par_r <- 0.2
# función a integrar
crec2 <- function(t, N, par_r){
        with(
            as.list(c(N, par_r)),{
                dx <- par_r*N
                return(list(dx))
            }
            )
}
# integración
out <- ode(y = N, times = time, func = crec2, parms = par_r)
# cambiando el formato del archivo de salida
ninteg <- as.data.frame(out)
colnames(ninteg) <- c("time","N")
# valores de N
head(ninteg)
##   time        N
## 1  0.0 100.0000
## 2  0.1 102.0202
## 3  0.2 104.0812
## 4  0.3 106.1838
## 5  0.4 108.3288
## 6  0.5 110.5172
# gráfica
ggplot(data = ninteg, aes(x=time, y=N)) + geom_point()

Ajuste del modelo a la realidad

Sabemos que las poblaciones (excepto la humana por un buen tiempo) no crecen de manera exponencial todo el tiempo. Vamos a introducir un elemento de control en la ecuación anterior, que reduce el crecimiento a medida que la población se acerca a cierto valor (capacidad de acarreo, \(K\)).

La nueva ecuación diferencial sería: \[\frac{dN}{dt} = rN(1 - \frac{N}{K})\]

La solución analítica de esta ecuación diferencial es: \[N(t) = \frac{KN_0}{N_0 + (K - N_0)e^{-rt}}\]

Como ven, obtener la solución a la ecuación diferencial se va complicando, y en muchos casos no hay una solución analítica y solo aproximaciones numéricas. Por esto para la resolución de las ecuaciones y los posteriores análisis del modelo y simulaciones, es mejor utilizar directamente deSolve u otro similar.

# tiempo
time <- seq(from=0, to=30, by = 0.1)
# valor inicial de N
N <- 100
# valor del parámetro a, velocidad de crecimiento
parametros <- c(r = 0.4, K = 1000)
# función a integrar
crec3 <- function(t, N, parametros){
        with(as.list(c(N, parametros)),{
                dx <- r*N*(1 - N/K)
                return(list(dx))
            }
            )
}
# integración
out <- ode(y = N, times = time, func = crec3, parms = parametros)
# cambiando el formato del archivo de salida
ninteg <- as.data.frame(out)
colnames(ninteg) <- c("time","N")
# valores de N
head(ninteg)
##   time        N
## 1  0.0 100.0000
## 2  0.1 103.6581
## 3  0.2 107.4340
## 4  0.3 111.3303
## 5  0.4 115.3497
## 6  0.5 119.4947
# gráfica
ggplot(data = ninteg, aes(x=time, y=N)) + geom_point()

Tasas de cambio instantáneas en un sistema de reacción tardía

El modelo logístico, aunque más realista en su descripción mecanicista del funcionamiento de las poblaciones, no refleja el comportamiento oscilante que muchas veces se observa en organismos superiores. En el modelo logístico básico, la respuesta al aumento de la población mediante el ajuste con \(K\), es instantánea, lo cual en los sistemas biológicos es muy poco probable. Para esto, se puede incluir una función de retraso, que realice la integración utilizando valores previos de \(N\). Utilizaremos la función lagvalue y el integrador dede para un modelo con retraso (delay).

# tiempo
time <- seq(from=0, to=30, by = 0.1)
# valor inicial de N
N <- 10
# valor del parámetro a, velocidad de crecimiento
parametros <- c(r = 0.4, K = 1000)
# función a integrar
crec3 <- function(t, N, param){
        with(as.list(c(N, parametros)),{
  tlag <- t - 2.3
  if (tlag < 0)
    Nlag <- N
  else 
    Nlag <- lagvalue(tlag)
                dx <- r*N*(1 - Nlag/K)
                return(list(dx))
            }
            )
}
# integración
out <- dede(y = N, times = time, func = crec3, parms = NULL)
# cambiando el formato del archivo de salida
ninteg <- as.data.frame(out)
colnames(ninteg) <- c("time","N")
# valores de N
head(ninteg)
##   time        N
## 1  0.0 10.00000
## 2  0.1 10.40387
## 3  0.2 10.82386
## 4  0.3 11.26062
## 5  0.4 11.71479
## 6  0.5 12.18705
# gráfica
ggplot(data = ninteg, aes(x=time, y=N)) + geom_point()