Año 2 · Módulo 13

127.1. Fluctuaciones demográficas en poblaciones pequeñas

En una población numerosa, la trayectoria media puede sugerir una dinámica regular. En una población pequeña, una sucesión desfavorable de muertes puede provocar extinción incluso cuando la tasa individual de reproducción supera a la de mortalidad. Los procesos de nacimiento y muerte hacen calculable esa diferencia. Su estado es el número entero de unidades presentes y sus mecanismos se expresan mediante tasas de transición que dependen del tamaño poblacional.

Las fluctuaciones son especialmente relevantes cuando los conteos son pequeños. Una muerte es una variación relativa del 12,5 % si quedan ocho individuos y del 50 % si quedan dos. Una variable continua no experimenta, por definición, saltos indivisibles de una unidad. Tampoco reproduce por sí sola el carácter absorbente del cero; una vez extinguida una población cerrada, sin inmigración ni mecanismos de aparición espontánea, el modelo no dispone de individuos capaces de reproducirse. En una aproximación determinista, el tamaño puede acercarse a cero sin alcanzarlo nunca. Ambas descripciones pueden ser útiles, siempre que indiquemos sus respectivos dominios.

El problema no se limita a la ecología. El número de moléculas de una especie química, los portadores de una infección en sus etapas iniciales, las células de un clon raro y los defectos en un material pueden formularse mediante conteos enteros. La pregunta puede dirigirse al número promedio de unidades, a la posibilidad de extinción, al tiempo que transcurre hasta ella o a la distribución completa de tamaños. Para estas últimas preguntas necesitaremos conservar la aleatoriedad de los acontecimientos individuales.

Para una población de tamaño nn, una tasa total de nacimientos proporcional a nn describe una situación distinta de una tasa que disminuye a causa de la competencia por recursos. En el primer caso pueden obtenerse fórmulas cerradas para el crecimiento medio y la probabilidad de extinción; en el segundo, la dependencia de las tasas respecto de nn modifica esas probabilidades. La interpretación requiere mantener separados el modelo lineal, las restricciones de densidad y la eventual introducción de una capacidad finita.

La ecuación maestra debe conservar el carácter absorbente del cero en un sistema cerrado, además de normalizar las probabilidades y mantener no negativas todas las tasas. Estas propiedades son controles independientes para detectar errores en soluciones numéricas y simulaciones de eventos.

127.2. El proceso de nacimiento y muerte como cadena de Markov

Sea N(t)N(t) el número de individuos presentes en el instante tt. Su espacio de estados es S={0,1,2,…}\mathcal S=\{0,1,2,\ldots\}. Desde un estado nn solo permitimos dos saltos elementales, el nacimiento n→n+1n\to n+1 y la muerte n→n−1n\to n-1. Denotaremos por bnb_n y dnd_n sus respectivas tasas. Para n=0n=0, la muerte no es posible, de modo que d0=0d_0=0. Si estudiamos una población cerrada sin inmigración, también imponemos b0=0b_0=0; entonces el estado cero es absorbente.

La interpretación infinitesimal requiere formular con cuidado las probabilidades. Para un intervalo pequeño h>0h>0 y condicionado a N(t)=nN(t)=n, tenemos

P(N(t+h)=n+1∣N(t)=n)=bnh+o(h),P(N(t+h)=n−1∣N(t)=n)=dnh+o(h),P(N(t+h)=n∣N(t)=n)=1−(bn+dn)h+o(h).\begin{aligned} P(N(t+h)=n+1\mid N(t)=n)&=b_nh+o(h),\\ P(N(t+h)=n-1\mid N(t)=n)&=d_nh+o(h),\\ P(N(t+h)=n\mid N(t)=n)&=1-(b_n+d_n)h+o(h). \end{aligned}

Según estas expresiones, la probabilidad de un salto determinado es, en primera aproximación, proporcional a hh; la tasa no fija el número exacto de nacimientos por unidad de tiempo. En un intervalo finito pueden ocurrir varios eventos, siempre que las tasas sean adecuadas para construir un proceso no explosivo. La tasa total de salida del estado nn es an=bn+dna_n=b_n+d_n y la permanencia allí, si an>0a_n>0, tiene distribución exponencial de parámetro ana_n.

Con la convención de generadores por filas, los elementos no nulos del generador QQ satisfacen

qn,n+1=bn,qn,n−1=dn,qn,n=−(bn+dn).q_{n,n+1}=b_n,\qquad q_{n,n-1}=d_n,\qquad q_{n,n}=-(b_n+d_n).

Es una matriz tridiagonal, infinita en el modelo sin límite superior. Por su estructura, cada salto individual cambia el tamaño en una sola unidad. Si una observación muestra variaciones mayores durante una unidad de tiempo, lo probable es que hayan ocurrido varios saltos entre dos mediciones. La discretización observacional de la trayectoria no coincide con la discreción del mecanismo.

Sea pn(t)=P{N(t)=n}p_n(t)=P\{N(t)=n\}. En ausencia de flujos externos, la ecuación maestra se escribe, para n≥1n\geq1,

dpndt=bn−1pn−1+dn+1pn+1−(bn+dn)pn.\frac{dp_n}{dt}=b_{n-1}p_{n-1}+d_{n+1}p_{n+1}-(b_n+d_n)p_n.

La ecuación para el estado cero, cuando b0=0b_0=0, es

dp0dt=d1p1.\frac{dp_0}{dt}=d_1p_1.

Por ello, la probabilidad acumulada de extinción p0(t)p_0(t) no disminuye; la masa probabilística que entra en cero no regresa a estados positivos. La evolución completa es un balance entre flujos de llegada y salida de cada tamaño poblacional. Al sumar las ecuaciones desaparecen los flujos internos y se conserva la probabilidad total, siempre bajo las condiciones regulares necesarias para un espacio infinito de estados.

127.3. Nacimiento puro, crecimiento medio y dispersión de trayectorias

Consideremos primero un modelo en el que cada individuo da origen a otro con tasa constante λ>0\lambda>0, sin mortalidad. Si hay nn individuos, hay nn oportunidades de reproducción. La tasa total de nacimiento es, por tanto,

bn=λn,dn=0.b_n=\lambda n,\qquad d_n=0.

Este proceso, conocido como proceso de Yule, parte por ejemplo de N(0)=1N(0)=1. Cada nacimiento incrementa el número de futuros reproductores. El tiempo de espera hasta el siguiente salto es exponencial con tasa λn\lambda n; a medida que la población crece, los intervalos tienden a acortarse. A diferencia del proceso Poisson homogéneo de C125, no hay un reloj con tasa constante independiente del estado.

Para una sola unidad inicial, la distribución exacta es geométrica sobre {1,2,…}\{1,2,\ldots\}, con

P(N(t)=n)=e−λt(1−e−λt)n−1,n≥1.P(N(t)=n)=e^{-\lambda t}\bigl(1-e^{-\lambda t}\bigr)^{n-1},\qquad n\geq1.

Podemos comprobar que sus probabilidades suman uno mediante la serie geométrica, con razón 1−e−λt1-e^{-\lambda t}. Su esperanza y varianza resultan

E[N(t)]=eλt,Var⁡[N(t)]=eλt(eλt−1).E[N(t)]=e^{\lambda t},\qquad \operatorname{Var}[N(t)]=e^{\lambda t}\left(e^{\lambda t}-1\right).

La media coincide con la solución de dn/dt=λndn/dt=\lambda n con condición inicial uno. Las realizaciones, sin embargo, no son curvas exponenciales suaves; cada población permanece constante durante intervalos aleatorios y crece mediante saltos. Algunas realizaciones acumulan reproducciones tempranas; otras siguen con un solo individuo durante un tiempo considerable. La variabilidad se amplifica con el crecimiento y la desviación típica llega a ser comparable con la media.

Si comenzamos con mm individuos independientes e idénticos, cada uno origina su propio árbol de descendientes. El total es la suma de mm procesos de nacimiento puro y

E[N(t)]=meλt,Var⁡[N(t)]=meλt(eλt−1).E[N(t)]=m e^{\lambda t},\qquad \operatorname{Var}[N(t)]=m e^{\lambda t}(e^{\lambda t}-1).

La variación relativa típica disminuye inicialmente al aumentar el número de fundadores. Es un primer indicio de por qué los modelos macroscópicos pueden funcionar mejor en poblaciones extensas, donde las sumas de múltiples contribuciones independientes tienden a presentar fluctuaciones relativas menores. Esa afirmación supone independencia y una escala adecuada; si hay correlaciones o sucesos compartidos, no puede darse por supuesta.

En el modelo de Yule no puede producirse extinción porque no existe mecanismo de muerte. Tampoco ocurre una explosión de infinitos nacimientos en tiempo finito, porque la suma de los tiempos medios sucesivos, ∑n=1∞1/(λn)\sum_{n=1}^{\infty}1/(\lambda n), diverge. Este control es necesario cuando las tasas crecen con el tamaño, pues una regla que incremente la intensidad demasiado deprisa puede producir infinitos eventos acumulados en un tiempo finito.

Doce trayectorias escalonadas distintas de una población de Yule crecen alrededor de una curva media exponencial. A la derecha la distribución de tamaños es geométrica.
Figura 127.1. Doce trayectorias de nacimientos puros (semilla 12702) y la media exacta e^{0,7t}. A la derecha, ley geométrica en t = 2,2; se indica explícitamente la probabilidad de la cola no dibujada (n ≥ 18). Ampliar figura ↗

127.4. Nacimiento y muerte lineales e insuficiencia del signo del crecimiento

Incorporemos ahora mortalidad individual con tasa μ>0\mu>0. El proceso lineal de nacimiento y muerte queda definido por

bn=λn,dn=μn,n≥1.b_n=\lambda n,\qquad d_n=\mu n,\qquad n\geq1.

El estado cero continúa siendo absorbente. Cada individuo puede reproducirse o morir y una trayectoria termina cuando ya no queda ninguno. Antes del siguiente salto, el tiempo de permanencia en nn sigue una exponencial de tasa (λ+μ)n(\lambda+\mu)n. Condicionado a que ocurra un evento, la probabilidad de nacimiento es λ/(λ+μ)\lambda/(\lambda+\mu) y la de muerte es μ/(λ+μ)\mu/(\lambda+\mu). Estas probabilidades son independientes de nn solo por la elección particular de tasas lineales; en los modelos regulados por densidad dependen de nn.

El comportamiento promedio parece sencillo. Usando el generador aplicado a la función f(n)=nf(n)=n, obtenemos

ddtE[N(t)]=E[bN(t)−dN(t)]=(λ−μ)E[N(t)].\frac{d}{dt}E[N(t)]=E[b_{N(t)}-d_{N(t)}] =(\lambda-\mu)E[N(t)].

Si N(0)=mN(0)=m, la solución es

E[N(t)]=me(λ−μ)t.\boxed{E[N(t)]=m e^{(\lambda-\mu)t}.}

El signo de r=λ−μr=\lambda-\mu distingue crecimiento medio, balance crítico y decrecimiento medio. Sin embargo, el resultado es una esperanza sobre muchas realizaciones, incluidas las ya extinguidas. Si λ>μ\lambda>\mu, la media aumenta aunque una fracción de las poblaciones termine para siempre en cero; las trayectorias supervivientes pueden crecer lo suficiente para elevar el promedio, mientras otras contribuyen exactamente cero. El crecimiento de la media no equivale a la supervivencia de una población concreta.

La ecuación para la varianza también puede obtenerse del generador, ahora aplicado a n2n^2. Para r≠0r\neq0, partiendo de mm individuos, resulta

Var⁡[N(t)]=mλ+μrert(ert−1).\operatorname{Var}[N(t)] =m\frac{\lambda+\mu}{r}e^{rt}\left(e^{rt}-1\right).

Esta expresión permanece no negativa tanto si r>0r>0 como si r<0r<0, pues el factor (ert−1)/r(e^{rt}-1)/r es positivo. En el caso crítico λ=μ\lambda=\mu, el límite es

Var⁡[N(t)]=2mλt.\operatorname{Var}[N(t)]=2m\lambda t.

Una población crítica conserva media constante mm, pero su varianza crece linealmente en el tiempo. Cada trayectoria sigue experimentando nacimientos y muertes; parte de la masa se acumula en cero y otra parte se desplaza hacia tamaños mayores. La igualdad de tasas no describe un equilibrio estacionario de la distribución poblacional condicionada a estar viva.

Este modelo permite separar tres escalas. La tasa de un individuo describe el mecanismo local; la esperanza es el promedio de realizaciones; la probabilidad de extinción expresa un resultado global irreversible. Un modelo que solo informe rr omite las dos últimas piezas.

Curvas de probabilidad de extinción eventual según cociente de tasas y número de fundadores uno, dos y cinco; las tres alcanzan uno cuando mortalidad supera nacimiento.
Figura 127.2. Extinción eventual del proceso lineal de nacimiento y muerte sin límite de tamaño. Para μ < λ, q = μ/λ y para m fundadores independientes la probabilidad es q^m; con μ ≥ λ, la extinción ocurre casi seguramente. Ampliar figura ↗

127.5. Probabilidad y tiempos de extinción

Denotemos por qq la probabilidad de extinción eventual cuando comenzamos con exactamente un individuo. Suponemos λ>0\lambda>0, μ>0\mu>0 y que los descendientes evolucionan de manera independiente bajo las mismas tasas. El primer acontecimiento es muerte con probabilidad μ/(λ+μ)\mu/(\lambda+\mu) o nacimiento con probabilidad λ/(λ+μ)\lambda/(\lambda+\mu). Si ocurre muerte, la población se extingue de inmediato. Si ocurre nacimiento, quedan dos individuos y cada linaje debe extinguirse para que desaparezca toda la población. Por independencia, la probabilidad conjunta es q2q^2.

La condición de primer paso queda

q=μλ+μ+λλ+μq2.q=\frac{\mu}{\lambda+\mu} +\frac{\lambda}{\lambda+\mu}q^2.

Al multiplicar por λ+μ\lambda+\mu y reorganizar obtenemos λq2−(λ+μ)q+μ=0,\lambda q^2-(\lambda+\mu)q+\mu=0, o bien (q−1)(λq−μ)=0.(q-1)(\lambda q-\mu)=0.

Las raíces son q=1q=1 y q=μ/λq=\mu/\lambda. Si λ≤μ\lambda\leq\mu, la única raíz admisible en [0,1][0,1], o la raíz mínima en el caso crítico, es uno. Si λ>μ\lambda>\mu, existen dos soluciones matemáticas admisibles y la teoría de procesos de ramificación selecciona la menor solución no negativa de la ecuación de punto fijo. En consecuencia,

q={1,λ≤μ,μ/λ,λ>μ.\boxed{q=\begin{cases} 1,&\lambda\leq\mu,\\ \mu/\lambda,&\lambda>\mu. \end{cases}}

La ecuación algebraica admite ambas raíces, de modo que la elección de la menor requiere justificación. Una forma consiste en aproximar la probabilidad de extinción mediante probabilidades de extinción hasta una generación o un número finito de acontecimientos; la iteración desde cero converge a la menor solución. También puede demostrarse con métodos de martingalas o mediante la función generadora del proceso.

Si partimos de mm fundadores con linajes independientes, la probabilidad de extinción eventual es qmq^m. Tomemos λ=0,6\lambda=0{,}6 y μ=0,4\mu=0{,}4 por día. El crecimiento medio es positivo, r=0,2r=0{,}2 por día, pero desde un fundador la extinción eventual tiene probabilidad 2/32/3. Desde tres fundadores independientes, es (2/3)3=8/27(2/3)^3=8/27, aproximadamente 0,29630{,}2963.

La palabra «eventual» merece precisión. Cuando λ>μ\lambda>\mu, 1−q1-q es la probabilidad de supervivencia indefinida en un modelo ideal con recursos ilimitados y una población potencialmente sin cota superior; un ecosistema real no cumple esos supuestos. En un espacio finito de individuos y con absorción accesible desde cualquier estado, la probabilidad de extinción final puede ser uno incluso si existe crecimiento medio favorable a tamaños bajos.

La probabilidad de extinción eventual no informa cuándo ocurre. Para distinguir riesgo temprano de supervivencia prolongada definimos el tiempo de extinción

T0=inf⁡{t≥0:N(t)=0}.T_0=\inf\{t\geq0:N(t)=0\}.

Si la población comienza con un individuo, qt=P(T0≤t)q_t=P(T_0\leq t). En un proceso lineal con tasas λ\lambda y μ\mu, su expresión exacta, para λ≠μ\lambda\ne\mu, es

qt=μ(1−e−(λ−μ)t)λ−μe−(λ−μ)t.q_t=\frac{\mu(1-e^{-(\lambda-\mu)t})} {\lambda-\mu e^{-(\lambda-\mu)t}}.

Si λ=μ=γ>0\lambda=\mu=\gamma>0, la expresión límite es

qt=γt1+γt.q_t=\frac{\gamma t}{1+\gamma t}.

Para mm linajes independientes, la probabilidad de extinción de todos antes del instante tt es qtmq_t^m. Estas fórmulas suponen población cerrada, individuos equivalentes y ausencia de variaciones temporales en las tasas. Cuando cambia el ambiente o existe dependencia entre linajes, la factorización requiere revisión.

La fórmula crítica tiene una propiedad poco intuitiva. Aunque la extinción eventual es segura, la probabilidad de supervivencia hasta tt es

P(T0>t)=11+γt.P(T_0>t)=\frac{1}{1+\gamma t}.

Su descenso es algebraico. La esperanza del tiempo de extinción es entonces infinita, porque E[T0]=∫0∞P(T0>t) dtE[T_0]=\int_0^\infty P(T_0>t)\,dt diverge logarítmicamente. Una realización típica se extingue en un tiempo finito, pero la cola de los tiempos de supervivencia es tan extensa que el promedio matemático no es finito. Aquí sería un error describir el tiempo de extinción solo mediante su media.

En el caso supercrítico existe una probabilidad positiva de T0=∞T_0=\infty, de modo que su esperanza no condicionada también es infinita. El fenómeno es distinto del crítico, porque una parte de la masa corresponde a trayectorias que nunca se extinguen en el modelo ideal. Para describir los tiempos observables de las realizaciones que sí se extinguen debemos condicionar al evento T0<∞T_0<\infty, con el correspondiente cambio de distribución.

En el caso subcrítico, λ<μ\lambda<\mu, la extinción eventual tiene probabilidad uno. Para un fundador, la supervivencia decae asintóticamente de modo exponencial, con escala determinada por μ−λ\mu-\lambda. Cerca de la criticidad, la diferencia entre las tasas puede ser pequeña y los tiempos característicos muy extensos, aunque cada evento individual siga ocurriendo con tasas apreciables.

Conviene distinguir el tiempo de la última observación, el tiempo de la última muerte y el instante exacto de llegada a cero. En registros discretos puede conocerse solamente un intervalo donde ocurrió la extinción. Esta censura por intervalo altera la inferencia de T0T_0. Que un modelo tenga expresión analítica no exime de examinar cómo se obtuvieron los datos que lo alimentan.

Tres curvas de probabilidad acumulada de extinción: la supercrítica se aproxima a cuatro novenos, mientras que las críticas y subcríticas se aproximan a uno con velocidades diferentes.
Figura 127.3. Distribución acumulada del tiempo de extinción para un único fundador, en regímenes supercrítico, crítico y subcrítico. En el supercrítico, la asíntota es μ/λ = 4/9; las otras curvas tienden a uno. Ampliar figura ↗

127.6. Regulación demográfica y cuasiestacionariedad

El crecimiento indefinido del proceso lineal es una idealización. Si la reproducción depende de recursos y competencia, las tasas individuales pueden cambiar con el tamaño. Un modelo estocástico de tipo logístico conserva nacimientos lineales bn=λnb_n=\lambda n, pero asigna muertes

dn=μn+cKn(n−1),c>0,K>0.d_n=\mu n+\frac{c}{K}n(n-1),\qquad c>0,\quad K>0.

El primer término describe mortalidad basal; el segundo añade competencia entre pares, proporcional a n(n−1)n(n-1). Elegimos n(n−1)n(n-1), en lugar de n2n^2, porque un solo individuo no compite consigo mismo en este mecanismo específico. El parámetro KK introduce una escala de tamaño y cc determina la intensidad de la interacción competitiva. Sin calibración adicional, ninguno de los dos equivale a una capacidad de carga empírica.

La ecuación exacta de la media es

ddtE[N]=(λ−μ)E[N]−cKE[N(N−1)].\frac{d}{dt}E[N] =(\lambda-\mu)E[N]-\frac{c}{K}E[N(N-1)].

El término E[N(N−1)]E[N(N-1)] depende del segundo momento, pues E[N(N−1)]=E[N]2−E[N]+Var⁡(N)E[N(N-1)]=E[N]^2-E[N]+\operatorname{Var}(N). Si reemplazamos el término esperado por una función del promedio obtenemos una aproximación de cierre. Para grandes tamaños y fluctuaciones relativamente pequeñas, podemos aproximar E[N(N−1)]E[N(N-1)] mediante E[N]2E[N]^2, con las correcciones que correspondan, y llegar a una ecuación semejante a la logística.

La aproximación determinista resultante es

dndt≈(λ−μ)n−cKn2.\frac{dn}{dt}\approx(\lambda-\mu)n-\frac{c}{K}n^2.

Cuando λ>μ\lambda>\mu, esta ecuación tiene un equilibrio positivo cercano a K(λ−μ)/cK(\lambda-\mu)/c. El proceso estocástico, en cambio, mantiene accesible el estado cero. Puede fluctuar alrededor de un tamaño aparentemente estable durante mucho tiempo y finalmente extinguirse. Llamar «equilibrio» a la meseta observada no demuestra que exista una distribución estacionaria no trivial sobre poblaciones vivas. Una distribución estacionaria genuina puede concentrarse por completo en cero aunque las realizaciones supervivientes parezcan estabilizadas durante horizontes prácticos.

Este ejemplo también muestra el peligro de cerrar ecuaciones de momentos sin evaluar su efecto. Si la varianza es elevada, el término competitivo esperado será distinto del obtenido insertando el tamaño medio en la tasa de muerte. Una ecuación ajustada sobre promedios puede fallar precisamente al estimar riesgos de extinción, porque esos riesgos dependen de la probabilidad que alcanza los estados pequeños.

La intensidad de las fluctuaciones puede cambiar con KK, con el régimen de parámetros y con la cercanía al estado absorbente. Una aproximación continua por difusión puede describir correctamente fluctuaciones moderadas alrededor de un tamaño grande, pero reproducir mal el evento raro de caer a cero. Antes de confiar en una estimación de extinción conviene comprobar si el mecanismo discreto, el tratamiento de la frontera y la escala de la muestra siguen siendo relevantes.

Cuando un sistema tiene absorción final pero puede permanecer mucho tiempo en estados positivos, resulta útil estudiar la distribución condicionada a no haber alcanzado cero. Definimos

p~n(t)=P(N(t)=n∣T0>t),n≥1.\widetilde p_n(t) =P(N(t)=n\mid T_0>t),\qquad n\geq1.

Una distribución cuasiestacionaria es una ley sobre estados positivos que permanece invariante al evolucionar y volver a condicionar por supervivencia. Su existencia y unicidad no se presuponen; dependen del proceso concreto. Si existe y se inicia en ella, la supervivencia suele presentar una ley exponencial en el tiempo con una tasa de pérdida asociada al régimen cuasiestacionario. Este concepto da una forma precisa de estudiar poblaciones que parecen estabilizadas durante mucho tiempo aun cuando el destino absorbente siga disponible.

La condición de supervivencia también introduce selección estadística. Supongamos que seguimos cien poblaciones pequeñas bajo el mismo modelo. Transcurridos varios días, las más vulnerables pueden haberse extinguido. Si calculamos el tamaño promedio solo entre las supervivientes, obtendremos una cantidad distinta de E[N(t)]E[N(t)] sobre la cohorte original. El promedio condicionado omite precisamente las observaciones iguales a cero y, por ello, puede aumentar incluso si la media no condicionada decrece. Para interpretar series de poblaciones cultivadas, colonias microbianas o focos infecciosos resulta indispensable declarar qué unidades permanecieron en el denominador.

Como el proceso lineal es un proceso de ramificación, parte de la intuición se traslada a epidemias incipientes. Si cada portador origina nuevas infecciones con tasa aproximada β\beta y deja de transmitir con tasa γ\gamma, el número inicial de portadores puede aproximarse por un nacimiento-muerte lineal. El cociente R=β/γR=\beta/\gamma, en esta idealización de contactos efectivos y susceptibles prácticamente constantes, controla si la esperanza crece o disminuye. Cuando R>1R>1, la probabilidad de extinción eventual del brote aproximado desde un portador es γ/β\gamma/\beta; desde mm portadores independientes es (γ/β)m(\gamma/\beta)^m.

Este cociente no coincide necesariamente con todas las formulaciones del R0R_0 epidemiológico, cuya definición depende de la estructura de contactos, las duraciones infecciosas y la heterogeneidad. Tampoco es válido proyectar el nacimiento-muerte lineal durante toda una epidemia cerrada, porque al disminuir los susceptibles la tasa efectiva de transmisión deja de ser constante.

Hay otra precaución al transferir el lenguaje a sistemas humanos. El número de personas con un evento clínico, una hospitalización o una conducta no sigue automáticamente un proceso de nacimiento y muerte. Deben identificarse transiciones individuales observables, unidades de conteo, mecanismos de aparición y desaparición y supuestos de independencia. El lenguaje de «contagio», «resiliencia» o «extinción» puede ser una metáfora útil, pero no justifica científicamente el generador elegido. La estructura probabilística debe derivarse del fenómeno y de la calidad de sus datos.

127.7. Tasas, extinción y simulación en un caso finito verificable

Para unir la teoría con un cálculo reproducible estudiaremos una colonia con capacidad máxima idealizada de tres individuos. El espacio de estados es {0,1,2,3}\{0,1,2,3\}. Cada individuo se reproduce con tasa λ=1,2\lambda=1{,}2 por día y muere con tasa μ=0,8\mu=0{,}8 por día. Cuando la colonia alcanza tres individuos se suprimen los nacimientos; las muertes continúan. Es una regla de bloqueo, distinta de la regulación cuadrática de la sección anterior. Partimos de N(0)=2N(0)=2.

Las tasas de nacimiento son bn=1,2nb_n=1{,}2n para n=1,2n=1,2 y b3=0b_3=0. Las tasas de muerte son dn=0,8nd_n=0{,}8n. El generador por filas queda

Q=(00000.8−2.01.2001.6−4.02.4002.4−2.4).Q=\begin{pmatrix} 0&0&0&0\\ 0.8&-2.0&1.2&0\\ 0&1.6&-4.0&2.4\\ 0&0&2.4&-2.4 \end{pmatrix}.

Cada fila suma cero y los coeficientes fuera de la diagonal son no negativos. El estado cero es absorbente y puede alcanzarse desde todos los estados positivos mediante una secuencia de muertes. Como el espacio es finito y no existen otras clases cerradas, la extinción eventual ocurre con probabilidad uno. Esto contrasta con el modelo lineal sin límite de tamaño, en el que λ>μ\lambda>\mu permitía supervivencia indefinida con probabilidad positiva.

Para calcular la distribución en un instante tt, usamos p(t)=p(0)eQtp(t)=p(0)e^{Qt}. El tiempo medio de extinción puede hallarse sin simular. Sea un=E[T0∣N(0)=n]u_n=E[T_0\mid N(0)=n], con u0=0u_0=0. Las ecuaciones de primer paso son

−1=bn(un+1−un)+dn(un−1−un),n=1,2,3,-1=b_n(u_{n+1}-u_n)+d_n(u_{n-1}-u_n),\quad n=1,2,3,

con b3=0b_3=0. Si QTQ_T es la submatriz del generador correspondiente a estados positivos, el sistema se escribe

(−QT)u=1.(-Q_T)u=\boldsymbol 1.

La ecuación maestra hace evolucionar probabilidades, mientras que el sistema lineal de primer paso calcula tiempos esperados de llegada; ambos usan las mismas tasas, con incógnitas y condiciones distintas.

El siguiente programa calcula probabilidades exactas de extinción en tiempos elegidos, resuelve las esperanzas y las contrasta con trayectorias generadas mediante relojes exponenciales. Requiere NumPy y SciPy. El número de simulaciones se mantiene moderado para que el ejemplo pueda ejecutarse en un ordenador personal.

import numpy as np
from scipy.linalg import expm

rng = np.random.default_rng(127)
Q = np.array([
    [0.0, 0.0, 0.0, 0.0],
    [0.8, -2.0, 1.2, 0.0],
    [0.0, 1.6, -4.0, 2.4],
    [0.0, 0.0, 2.4, -2.4],
])
assert np.allclose(Q.sum(axis=1), 0.0)

p0 = np.array([0.0, 0.0, 1.0, 0.0])
for t in (1.0, 4.0, 12.0):
    pt = p0 @ expm(Q * t)
    print(f't={t:4.1f}: P(extinta)={pt[0]:.5f}, '
          f'E[N]={pt @ np.arange(4):.5f}')

T = Q[1:, 1:]
u = np.linalg.solve(-T, np.ones(3))
print('Tiempo medio a cero desde n=1,2,3:', u)

def una_extincion(n0=2):
    t, n = 0.0, n0
    while n > 0:
        nacimientos = 1.2 * n if n < 3 else 0.0
        muertes = 0.8 * n
        total = nacimientos + muertes
        t += rng.exponential(1.0 / total)
        n += 1 if rng.random() < nacimientos / total else -1
    return t

muestras = np.array([una_extincion() for _ in range(30_000)])
print('Promedio simulado desde n=2:', muestras.mean())
print('Promedio exacto desde n=2:', u[1])

A partir de dos individuos, la solución matricial da una probabilidad exacta de extinción de aproximadamente 0,171610{,}17161 al primer día, 0,596350{,}59635 al cuarto y 0,941220{,}94122 al duodécimo. Las respectivas esperanzas de tamaño son 1,723961{,}72396, 0,834310{,}83431 y 0,121490{,}12149. Estos números corresponden a la cohorte completa, contando como tamaño cero cada colonia extinguida. La solución del sistema de primer paso da u1=3,125u_1=3{,}125, u2=4,375u_2=4{,}375 y u3≈4,79167u_3\approx4{,}79167 días. Es posible que la colonia pase parte de su historia en el estado tres y, aun así, su destino eventual siga siendo cero con probabilidad uno.

Hay una comprobación adicional que no requiere calcular la exponencial matricial. Desde dos individuos no podemos llegar a cero mediante un solo salto; deben ocurrir, como mínimo, dos muertes sucesivas. Para tiempos muy pequeños, la primera muerte tiene tasa 1,61{,}6 y la segunda, desde el estado uno, tasa 0,80{,}8. La probabilidad de extinción comienza entonces con el término (1,6)(0,8)t2/2=0,64t2(1{,}6)(0{,}8)t^2/2=0{,}64t^2, mientras que su derivada inicial es cero. Una función numérica que predijera una probabilidad de extinción proporcional a tt desde n=2n=2 contradiría la estructura misma del proceso. Este tipo de control local es útil para comprobar código sin depender exclusivamente de una comparación global.

La salida exacta debe satisfacer probabilidad total uno y límites entre cero y uno. El promedio simulado fluctúa alrededor del valor analítico; no es razonable exigir igualdad decimal entre una muestra finita y una esperanza exacta. Si duplicamos las repeticiones, esperamos reducir el error de muestreo, aunque no en una proporción determinada para cada ejecución. La semilla fija garantiza la reproducibilidad de la secuencia pseudoaleatoria con una implementación compatible; la muestra sigue sin ser una demostración.

El ejemplo también permite comprobar un sesgo de interpretación. Si seleccionamos solo las colonias que continúan vivas a los doce días, su tamaño promedio condicionado puede superar ampliamente el promedio no condicionado. La diferencia procede de excluir las realizaciones absorbidas. Para estimar un riesgo real deberíamos registrar el tiempo de observación, el número inicial de colonias y las posibles pérdidas de seguimiento.

Izquierda probabilidad acumulada de extinción exacta y puntos Monte Carlo, derecha promedios del tamaño para toda la cohorte y solo entre supervivientes.
Figura 127.4. Solución de exp(Qt) y 22 000 simulaciones de extinción, semilla 12705. Las barras muestran intervalos normales de Monte Carlo aproximados al 95 %. A la derecha, el tamaño medio entre supervivientes es distinto del promedio original, que incluye las colonias extinguidas. Ampliar figura ↗
Diagrama de estados cero a tres: tasas de nacimiento 1,2 y 2,4; tasas de muerte 0,8, 1,6 y 2,4. A la derecha matriz generadora cuatro por cuatro.
Figura 127.5. Tasas del ejemplo finito de C127.7 y matriz generadora Q en unidades de día⁻¹. La fila 0 es nula; cada otra fila suma cero. Los nacimientos se bloquean al alcanzar n = 3. Ampliar figura ↗

127.8. Aplicabilidad de los modelos de nacimiento y muerte

Los procesos de nacimiento y muerte muestran por qué la dinámica de una población no queda resumida por una tasa neta. Dos modelos con idéntico valor de λ−μ\lambda-\mu pueden tener diferentes tasas totales de recambio, distribuciones de tiempos de espera y probabilidades de extinción en horizontes finitos. En el modelo lineal, la probabilidad de extinción eventual depende del cociente μ/λ\mu/\lambda. Para estudiar el destino de una población pequeña debemos distinguir el crecimiento medio, la variabilidad, la estructura del estado absorbente y el horizonte temporal de la decisión.

También cuenta la estructura del espacio de estados. Un proceso con recursos ilimitados puede sobrevivir indefinidamente con probabilidad positiva cuando los nacimientos dominan a las muertes. Un proceso finito cuya única clase cerrada es el estado cero acaba extinguiéndose con probabilidad uno. Una colonia regulada por densidad puede mostrar una estabilidad aparente durante muchísimo tiempo, aunque al final se absorba. Indicar solo «crecimiento positivo» o «capacidad de carga» deja indistinguibles estas situaciones.

Al contrastar modelos con observaciones es necesario reconstruir el mecanismo de muestreo. Un registro agregado de tamaños poblacionales no identifica por separado las tasas de nacimiento y muerte; aunque su diferencia esté bien estimada, puede haber muchas combinaciones compatibles con ella. Registrar eventos individuales y sus tiempos da información distinta de la de una serie de censos. Si las unidades interactúan, si existen cambios ambientales compartidos o si la población se compone de clases heterogéneas, la hipótesis de individuos independientes también deberá revisarse.

La probabilidad de extinción depende de la relación entre las tasas, del tamaño inicial y de la regulación poblacional. En el proceso lineal supercrítico, una población finita tiene probabilidad positiva de extinguirse aunque su esperanza crezca exponencialmente.

Las simulaciones del ejemplo permiten comparar la distribución completa de tamaños con su media. Al condicionar en la supervivencia, las colonias que permanecen activas forman una subpoblación seleccionada cuyo tamaño promedio no es el del conjunto inicial; esa selección debe figurar explícitamente en cualquier informe de resultados.

Referencias

  • Allen, L. J. S. (2010). An Introduction to Stochastic Processes with Applications to Biology (2.ª ed.). CRC Press. Procesos de nacimiento y muerte, ramificación y extinción.
  • Bailey, N. T. J. (1964). The Elements of Stochastic Processes with Applications to the Natural Sciences. Wiley. Formulación clásica de procesos poblacionales.
  • Harris, T. E. (1963). The Theory of Branching Processes. Springer. Probabilidades de extinción y procesos de ramificación.
  • Norris, J. R. (1997). Markov Chains. Cambridge University Press. Generadores, absorción y cadenas en tiempo continuo.
  • van Kampen, N. G. (2007). Stochastic Processes in Physics and Chemistry (3.ª ed.). Elsevier. Ecuaciones maestras, tamaños finitos y fluctuaciones.