Crecimiento exponencial (discreto y continuo). Aplicación en R.

1. CRECIMIENTO EXPONENCIAL


Crecimiento poblacional discreto: el cambio en el tamaño poblacional N(t) en el tiempo N(t+1), depende de 4 factores: nacimientos, muertes, inmigración y emigración.

N(t+1)=N(t)+B-D+I-E y el cambio en el tamaño poblacional deltaN=N(t+1)-N(t)=B-D+I-E
Si la población es cerrada: N(t+1)=N(t)+B-D y deltaN=B-D
donde B es el número de nacimientos por individuo, D la probabilidad de que un individuo muera durante un intervalo de tiempo específico y rd=B-D es el factor de crecimiento discreto.

Asume: el crecimiento es suavizado, por lo que los nacimientos y muertes ocurren continuamente.

Si las tasas de nacimientos y muertes dependen del tamaño poblacional presente, y sea b y d las tasas instantáneas de nacimientos y muertes, respectivamente, tenemos:

B=b*N(t) y D=d*N(t) (*)

llamadas tasas relativas o tasas per cápita.

Así, obtenemos: N(t+1)=(1+b-d)*N(t) y deltaN=(b-d)*N(t)=r*N(t) (=B-D )

donde lambda=(1+b-d)=N(t+1)/N(t) es la tasa de cambio en un período finito de tiempo o tasa finita de cambio per cápita (también denotada por G) y r=(b-d) es la tasa instantánea de crecimiento per cápita o tasa intrínseca de crecimiento o parámetro Malthusiano.
Nota: lambda es adimencional mientras que las unidades de r son de individuos generados por cada individuo en cada unidad de tiempo (ind/(ind*tiempo)).

Por tanto, la versión discreta del modelo de crecimiento exponencial viene dada por la ecuación en diferencias discreta:

N(t+1)=lambda*N(t)

En este caso las tasas de natalidad y mortalidad sean constantes en el tiempo, por lo que la población crece o decrece cada año en una misma proporción.

Resolviendo la ecuación para cualquier plazo temporal t obtenemos:

N(t)=lambda^t*N(0)

Crecimiento poblacional continuo: cambio en el tamaño poblacional (dN) durante un intervalo de tiempo infenitesimal (dt).

dN/dt=B-D

A pesar que el crecimiento continuo no es posible, debido a que los individuos no se reproducen en valores fraccionarios, el crecimiento poblacional de organismos de vida larga puede ser aproximadamente continuo. Este tipo de crecimiento se describe con ecuaciones diferenciales continuas.
Nuevamente, si estas las tasas instantáneas son constantes en el tiempo, y utilizando (*), obtenemos: dN/dt=(b-d)*N. El valor de r determina si la población crece exponencialmente (r>0), permanece constante (r=0) o disminuye hacia la extinción (r<0): style="text-align: center;">
dN/dt=r*N

Si integramos esta ecuación obtenemos la siguiente solución (necesaria por ejemplo para realizar predicciones):

N(t)= exp(r*t)*N(0)

donde N(0) y N(t) corresponden a los tamaños poblacionales en el tiempo t=0 y t=t, respectivamente. Su transformación logarítmica convierte la curva de crecimiento exponencial en una línea recta, donde la pendiente corresponde a r.

Entonces, el modelo continuo es equivalente al modelo discreto en pasos de tiempo infenitesimales, con: r=ln(lambda) o lambda=exp(r). Por ello, la tasa r también se llama tasa logarítmica de cambio per cápita o adecuación biológica de un individuo en un intervalo de tiempo determinado (identificado también con la letra R).

La tasa r puede calcularse a partir de la densidad poblacional:

r=ln(N(t+1))-ln(N(t))=ln(N(t+1)/N(t))

De manera equivalente podemos encontrar en la literatura la siguiente expresión:

r=ln(N(t))-ln(N(t-1))=ln(N(t)/N(t-1)).

Supuestos del modelo denso-independiente:

• Sin estructura espacial, no consideramos ni I ni E, la población está cerrada.
• b y d constantes. Esto implica recursos ilimitados para el crecimiento poblacional, su potencial para el crecimiento exponencial. Explica la naturaleza mutiplicativa del crecimiento poblacional y una retroalimentación (feed-back) positiva.
• Sin estructura genética, no existe variación genética en la población (todos los individuos tienen las mismas tasas de natalidad y mortalidad) o si existe, esta permanece constante en el tiempo (r representa la media de la tasa instantánea de crecimiento para los diferentes genotipos de la población).
• Sin estructura por edades o tamaños. No existen diferencias en b y d entre individuos debido a su edad o tamaño corporal.
• Crecimiento continuo sin retrasos temporales. Los individuos nacen y mueren continuamente.

Componentes estocásticos

• Estocasticidad ambiental:
Los modelos determinísticos representan una visión idealizada sobre un mundo simple y ordenado. Sin embargo, el mundo real es complejo e incierto. Podemos calcular 2 números que nos ayudaran a estimar la incertidumbre: la media y la varianza de los datos.
¿cómo podemos incorporar la incertidumbre? Si la tasa de aumento cambia de manera impredecible en el tiempo, podemos incorporar la incertidumbre en el parámetro r. Esta variabilidad asociada con años buenos y malos para el crecimiento poblacional, se conoce como estocasticidad ambiental. Así, el tamaño poblacional medio para una población creciendo con estocasticidad ambiental, viene dada por:

Mean(N(t))=N(0)*exp(mean(r)*t)

que se diferencia del modelo determinístico en que se utiliza la media de r para predecir la media de N(t). La varianza del tamaño poblacional en el tiempo t viene dada por:

sigma_N(t)^2=N(0)^2*exp(2*mean(r)*t)*[ exp(sigma_r^2* t)-1]

Si la varianza de r es nula, entonces la varianza de N(t) también será nula y volvemos al modelo determinístico de crecimiento. Asimismo, la población se extinguirá si la varianza de r es mayor a dos veces la media de r:

sigma_r^2>2*mean(r)

En el modelo estocástico el tamaño poblacional medio aumenta exponencialmente como función de mean(r).


• Estocasticidad demográfica
La estocasticidad ambiental no es la única fuente de variabilidad poblacional. Aún cuando r es constante, los tamaños poblacionales pueden fluctuar debido a la aleatoriedad asociada a los nacimientos y muerte, generalmente importante en poblaciones pequeñas. La probabilidad de nacimientos o muertes dependerá de la magnitud relativa de b y d:

P(nacimientos)=b/(b+d) y P(muertes)=d/(b+d)

Esta estocasticidad demográfica es análoga a la deriva genética donde la frecuencia alélica varía aleatoriamente en pequeñas poblaciones.
La varianza del tamaño poblacional en el tiempo t viene dada por:

Sigma_N(t)^2=2*N(0)*b*t si b=d.
Sigma_N(t)^2=N(0)*(b+d)*exp(r*t)*[exp(r*t)-1]/r si b=!d.

La probabilidad de extinción depende del tamaño relativo de d y b, así como del tamaño poblacional inicial:

P(extinción)=(d/b)^(N(t))

Aplicación en R (librería "primer")

Crecimiento exponencial (geométrico o Maltusiano) discreto

### 1) CRECIMIENTO DENSO-INDEPENDIENTE GEOMÉTRICO DISCRETO
#gráfico del tamaño poblacional
N <- c(1, 3, 9, 27, 81); year <- 2001:2005; plot(year, N)
#tasa de crecimiento
rates = N[2:5]/N[1:4]
#proyección de la población al futuro Nt=(lambda^t)*No ### lambda es la tasa finita de crecimiento (es la tasa per cápita de crecimiento si la población crece geométricamente) ### el factor de crecimiento discreto rd es lambda=(1+rd)
N0 <- 1; lambda <- 2; time <- 0:10; Nt <- N0 * lambda^time

#efectos del tamaño inicial poblacional
N0 <- c(10, 20, 30); lambda <- 2; time <- 0:4
Nt.s <- sapply(N0, function(n) n * lambda^time) #calculamos el tamaño poblacional cada vez aplicando la función (n*lambda^time) a cada elemento de No
matplot(time, Nt.s, pch = 1:3) #notar que cambia el intercepto y pendiente (pero no cambia la pendiente en escala logarítmica: cambia la tasa absoluta de crecimiento (N2-N1) pero no la tasa relativa de crecimiento (N2/N1)
matplot(time, Nt.s, log = "y", pch = 1:3) #cuando tomamos logaritmos a ambos lados, las tasas (lambda) se vuelven diferencias, lo cual resulta en lineas paralelas en el gráfico

#efectos de diferentes tasa de crecimiento per cápita
N0 <- 100;  time <- 0:3;  lambdas <- c(0.5, 1, 1.5) #cuando lambda<1 la población se hunde y si lambda>1 la población crece. (It will always be true that when lambda > 1 and t > 1, lambda^t > lambda. It will also be true that when lambda< 1 and t > 1, lambda^t < lambda). cuando lambda=1 la población no cambia porque 1^t=1 para todo t.
N.all <- sapply(lambdas, function(x) N0 * x^time)
matplot(time, N.all, xlab = "Years", ylab = "N", pch = 1:3);  abline(h = N0, lty = 3);  text(0.5, 250, expression(lambda > 1), cex = 1.2);  text(0.5, 20, expression(lambda < 1), cex = 1.2)

#tasa de crecimiento medio (veremos cómo la media aritmética no es conveniente)
#sea N0=100; N1=N0*(0.5); N2=N1*(1.5), es decir N2=N0*R1*R2 donde Nt+1/Nt es una variable aleatoria llamada R. ¿cómo calcular una media para elementos que se mutliplican?    #Si Rmeant= R1*R2 entonces(Rmeant)^(1/t)=(R1*R2*...*Rt)^(1/t) entonces Rmean=(R1*R2*...*Rt)^(1/t) esto se llama la media geométrica:  Rmean=(prod(Ri))^(1/t) #media geométrica (las proyecciones basadas en la media geométrica R son menores que la media aritmética)   #comparamos las medias aritmética y geométrica t=5;  SS6=sparrows[1:(t + 1), ]
SSgr=SS6$Count[2:(t + 1)]/SS6$Count[1:t];  lam.A=sum(SSgr)/t; lam.G=prod(SSgr)^(1/t)
N0=SS6$Count[1];  plot(0:t, SS6$Count, ylab = "Projected Population Size"); lines(0:t, N0 * lam.A^(0:t), lty = 2); lines(0:t, N0 * lam.G^(0:t), lty = 1); legend(0, 70, c("Arithmetic Ave.", "Geometric Ave."),title = "Projections Based On:", lty = 2:1, bty = "n", xjust = 0)



Crecimiento exponencial (geométrico o Maltusiano) continuo

### 2) CRECIMIENTO DENSO-INDEPENDIENTE EXPONENCIAL CONTINUO    (ej: para especies asincrónicas como E.coli)
#crecimiento geométrico: cambio proporcional de la población en un intervalo de tiempo finito especificado
#crecimiento exponencial: cambio proporcional de la población instantáneo (en un instante !)
#sea una población E.coli con lambda=1.5, recordando que lambda=(1+rd) entonces N1=N0*(1+0.5) (rd=0.5) y N1/N0=(1+rd/n)^n. Cuando n tiende a infinito
#una aproximación numérica al límite de la expresión (1+rd/n)^n:
n <- 0:100; N0 <- 1; rd <- 1 #sea lambda=2 entonces rd=1
N1 <- N0 * (1 + rd/n)^n #calculamos la expresión para grandes valores de n (n= el número de divisiones dentro de un año), entonces rd/n es la tasa finita de crecimiento durante esos pasos de tiempo fraccional.
plot(n, N1/N0, type = "l"); text(50, 2, "For n = 100,"); text(50, 1.6, bquote((1 + frac("r"["d"], "n"))^"n" == .(round(N1[101]/N0, 3)))) #Indeed, if a population grew in a discrete annual step Nt+1 = Nt (1 + rd), the same rd, divided up into many small increments, would result in a much larger increase.
#el crecimiento tiene una expresión matemática simple: exponencial (e). Es decir: lim(1+r/n)^n=exp(r)
#cuando la población crece geométricamente con pasos de tiempo infinitesimales, decimso que la población crece exponencialmente y lo representamos como Nt=No*exp(rt) donde r es la tasa de crecimiento per cápita o la tasa intrínseca de cercimiento.

#proyección del tamaño poblacional con crecimiento exponencial continuo
r <- c(-0.03, -0.02, 0, 0.02, 0.03); N0 <- 2; t <- 1:100 #elegimos 5 valores diferentes para r
cont.mat <- sapply(r, function(ri) N0 * exp(ri * t))   #graficamos en escala aritmética y logarítmica
layout(matrix(1:2, nrow = 1)); matplot(t, cont.mat, type = "l", ylab = "N", col = 1); legend("topleft", paste(rev(r)), lty = 5:1, col = 1, bty = "n", title = "r"); matplot(t, cont.mat, type = "l", ylab = "N", log = "y", col = 1)

#derivando la derivada del tiempo. Podemos diferenciar la ecuación Nt=N0*exp(rt) respecto al tiempo para pbtener una ecuación diferencial para la tasa de crecimiento poblacional instantánea:  d(XY)/dt=dX/dt*Y+dY/dt*X, entonces dN0*exp(rt)/dt=dN0*exp(r)^t*N0 (con derivada de una constante es cero y derivada de a^t es ln(a)*a^t) y dN0*exp(rt)/dt=0*exp(r)^t+ln(exp(r))*exp(r)^t*N0  (con ln(e)=1 y N0*exp(rt)=N), entonces dN/dt=rN
#para describir la tasa de cambio poblacional podemos utilizar lambda (ej. 1.1 quiere decir que hay un 10% de incremento poblacional) o r (pero es adimensional) o el "doubling time": 2N0=N0*exp(rt) entonces t=ln(2)/r
m.time <- function(r, m = 2) { log(m)/r }; rs <- c(0, 1, 2); m.time(rs) #evaluamos diferentes valores para m (cualquier múltiplo arbitrario de N0) y r
#relacionar lambda con r: si asumimos un crecimiento exponencial constante y geométrico, podemos calcular r desde los datos. Sea Nt=N0*exp(rt)  entonces  ln(Nt)=ln(N0)+rt por lo que r es la pendiente de la relación lineal entre ln(Nt) y el tiempo t, y ln(N0) es el intercepto con el eje y. Entonces lambda=exp(r) o ln(lambda)=r.

####### Los modelos DI asumen ("solo" según Steven) que N crece con tasa de crecimiento per cápita constante sobre intervalo(s) de tiempo, es decir, no existe relación entre N y la tasa de crecimiento.
####### Para explicar fenómenos complejos (e.g. lambda) en término de sus determinantes (e.g. B y D) podemos ver los modelos como: Nt+1=Nt+B*Nt-D*Nt, es decir lambda=1+(B-D) y rd=B-D, con B el número de nacimientos por individuo y D la probabilidad de que un individuo muera durante el intervalo de tiempo especificado. Incluso r=b-d las tasas instantáneras per cápita.

####### Modelando con datos: dinámicas simuladas. Para predecir la población en el futuro podemos utilizar distintas aproximaciones: 1) predicción determinística (utiliza Rmean) o 2) simulación de la dinámica poblacional.
#(http://www.mbr-pwrc.usgs.gov/bbs/  We will: 1. look at and collecting the data (annual R’s),2. simulate one projection,3. scale up to multiple simulations,4. simplify simulations and perform 1000’s, and 5. analyze the output.
names(sparrows); attach(sparrows);plot(Count ~ Year, type = "b")
obs.R<-Count[-1]/Count[-length(Count)]; plot(obs.R ~ Year[-length(Count)]); abline(h = 1, lty = 3) #calculamos Rt=Nt+1/Nt
years<-50; set.seed(3); sim.Rs<-sample(x=obs.R,size=years,replace=TRUE) #haremos simulaciones asumiendo que los Rs son representativos de R en el futuro y cada uno es igualmente probable de ocurrir. Remuestreamos los Rs observados con reemplazamiento para cada año de la simulación, lo que determina una realización aleatoria de una trayectoria poblacional posible. Decidimos qué cantidad de año queremos simular (years).
output <- numeric(years + 1);output[1] <- Count[Year == max(Year)];for (t in 1:years) output[t + 1] <- { output[t] * sim.Rs[t] };plot(0:years, output, type = "l") #Nt+1=R*Nt. El gráfico representa un posible resultado de la trayectoria si asumimos que R es igualmente probable dado las obsrevaciones rt. como una única trayectoria nos dice poco, debemos replicar el proceso para examinar la distribución de los resultados (momentos, etc)
sims = 10; sim.RM <- matrix(sample(obs.R, sims * years, replace = TRUE), nrow = years, ncol = sims) #múltiples simulaciones (sims=número de simulaciones realizadas)
output[1] <- Count[Year == max(Year)]; outmat <- sapply(1:sims, function(i) { for (t in 1:years){output[t + 1] <- output[t] * sim.RM[t,i]}; output });matplot(0:years, outmat, type = "l", log = "y")
#muchas simulaciones con una función
PopSim <- function(Rs, N0, years = 50, sims = 10) { sim.RM = matrix(sample(Rs, size = sims * years, replace = TRUE), nrow = years, ncol = sims); output <- numeric(years + 1); output[1] <- N0; outmat <- sapply(1:sims, function(i) { for (t in 1:years){output[t + 1] <- round(output[t] * sim.RM[t, i], 0)}; output }); return(outmat) }
#análisis de los resultados
N.2053 <- output[51, ]; summary(N.2053, digits = 6); quantile(N.2053, prob = c(0.0275, 0.975)) #ej: la distribución es asimétrica a la derecha (larga cola a la derecha) como una distribución lognormal (común en datos ecológicos, por lo cual se dice que la unidad "natural" de la población es log(N) ma´s que N)
hist(N.2053, main = "N"); hist(log10(N.2053 + 1), main = "log(N+1)"); abline(v = log10(quantile(N.2053, prob = c(0.0275, 0.975))+ 1), lty = 3) #graficamos los datos en escala original y utilizando logaritmos más uno, las líneas indican los cuantiles 2.5 y 97.5%.
logOR <- log(obs.R); n <- length(logOR); t.quantiles <- qt(c(0.025, 0.975), df = n - 1); #comparando las simulaciones con las proyecciones determinísicas podemos encontrar la t-distribución al 95% basado en los intervalos de confianza para la media geométrica de R. Pero como el logaritmo de la media geométrica de R es la media aritmética de logR, podemos calcular los intervalos de confianza t-basados en logR, destransformandolos para proyectar la población con bandas sup e inf. Tomamos logartitmos, calculamos los grados de libertar y los cuantiles relevantes para la distribución t.
se <- sqrt(var(logOR)/n); CLs95 <- mean(logOR) + t.quantiles * se #calculamos el error estándar de la media y los intervalos de confianza al 95% para logR
R.limits <- exp(CLs95); R.limits #destransformamos para tener R
N.Final.95 <- Count[Year == max(Year)] * R.limits^50; round(N.Final.95) #hacemos la proyección y vemos que el límite inferior para la proyección determinística es la misma (extinción) que para la simulación, pero el límite superior difiere. debemos estudiar el supuesto (asumimos que log R puede aproximarse con una t-distribución).
qqplot(qt(ppoints(n), df = n - 1), scale(logOR)); qqline(scale(logOR))        #comparamos los logR de los valores teóricos con la t-distribución. Scalamos logOR para hacer las comparaciones más claras. Se observa qu elos datos están más centrados de lo esperado y que existen 3 outliers

Referencias

• Berryman, A. A. 1999. Principles of population dynamics and their application. Stanley Thornes, UK.
• Gotelli, N.J. 2008. A Primer of Ecology. 4th edition. Sinauer Associates, Inc., Sunderland, MA.
• Stevens. M.H.H. 2009 A Primer of Ecology with R. Use R! Series. Springer.