General

Gompertz y Kermack-McKendrick

La epidemiologia de salón durante la pandemia se ha dividido en dos bandos: los que gustaban de modelar los casos diarios usando un modelo basado en dinamica de compartimentos SIR, SEIR, etc y los que preferian las funciones genericas de evolucion temporal como la ecuacion Logistica o la de Gompertz.

Desde el punto de vista de resolucion de ecuaciones diferenciales, el modelo de dos compartimentos S-I y la evolucion por funcion logistica son la misma cosa. La funcion que nos da el numero total de infectados sigue la ecuacion \(y’ = \beta y (N-y)\), donde los ceros del lado derecho nos indican entre que valores se mueve la funcion: desde un inicial de cero infectados hasta el total de la poblacion.

La funcion de Gompertz retoca el termino “de finalizacion” de la evolucion temporal, $N-y$, suavizando su llegada a cero con un comportamiento logaritmico

$$ y’ = -b y \log \frac y k$$

usando el “final size” $k$ que puede o no coincidir con el total de poblacion $N$.

El modelo SIR de Kermack-McKendrick, en su version mas simple, no se suele presentar con el “compartimento” de infectados totales (que viene a ser I+R), pero si lo usamos vemos que la sustitucion es en el factor que Gompertz deja quieto: se añade un comportamiento logaritmico al primer factor,

$$y’ = (\beta y + \alpha \log (1-\frac y N)) (N-y)$$

lo que crea un nuevo cero $k$ entre el trivial y el total; la solucion de la “final size relation”

$$ – \frac \alpha \beta \log (1-\frac k N) = k $$

de manera que ahora el tamaño de la epidemia evoluciona desde cero hasta $k$, asumiendo por supuesto que la solucion de esta ecuacion es real y menor que el total de poblacion $N$.

En resumen: tanto Gompertz como SIR introducen un término logarítmico en el cero de finalización, pero lo hacen de diferente manera. He preguntado en Math.SE si existe alguna forma conocida de transformar uno en otro, pero no creo que me den ninguna respuesta. Tanto en los ceros como en el maximo las dos expresiones se comportan de manera muy distinta.

de Vries con más digitos

GK Ottarsson publicó hace un par de años una breve nota sobre soluciones de ecuaciones de recursion en la que daba unos cuantos digitos mas para la alpha calculada por Hans de Vries, quedando

137.035999095829700489647400

que es un poco mas de los 1/0.00729735256865385342269 = 137.035999095829700489 que calculaba la nota original.

Recordemos que el original comparaba contra CODATA 2004. No se ha añadido mucha más precision desde entonces, y la cantidad ha envejecido bien, si consideramos que CODATA 2019 tiene un valor 137.035999083(20). Aunque la tension de los ultimos articulos apuntaba en la direccion contraria.

El factor de crecimiento siempre tiende a uno en una función polinómica

A ver, si tienes una función que crece como \(f(x) = b x^n\), y te pones a calcular el factor de incremento \(F(x) = f(x+a)/f(x)\), es obvio que tiende a uno:

$$F(x \to \infty) = (1+a/x)^n$$

Por supuesto, cualquier regimen exponencial \(f(x) = b e^{k x^n}\) tiene un factor de incremento \(F(x) = e^{k [(x+a)^n – x^n]}\) que se va a infinito para \(n>1\) mientras que para \(n=1\) tiende a $$F(x \to \infty) =e^{k a}$$ un límite que, al ser el exponente positivo, es siempre mayor que uno.

Así, la convergencia a un factor de crecimiento unidad tan solo nos indica que hemos abandonado el régimen exponencial. La “fase de ralentización”. De hecho, ni siquiera tenemos garantizado haber entrado en el regimen polinomico, podriamos estar en el subexponencial, n <1 del caso anterior. Por ejemplo exp(sqrt(x)).

Extensiones del modelito no compartimental de ayer

Primero un par de avisos de interes general:

Mirando esos apuntes y modelos se ve que los modelos compartimentales tienen muchas posibilidades, especialmente cuando comienzas a permitir que una valvula tenga dependencia funcional arbitraria, y no simplemente productos.

Pero bueno, ya que estaba con la idea no compartimental de ayer, voy a dejar aun anotadas dos o tres cosas que se me quedaron en el cuaderno.

La idea de poner un tiempo de incubacion gaussiano o con concentracion de delta de dirac no parece contraria a la intuicion, pero hacer eso con el tiempo de resolucion/removal no parece tan logico, habria que ir poniendo otras distribuciones d(t), Weibulls, cosas asi. Ahora, a falta de otro criterio conviene que podamos hacer analiticamente la integral, asi que vamos a hacer algo que suena cutre pero no deja de ser parecido a lo de los compartimentos: asumir que la probabilidad de detectar al infectado funciona con una distribucion exponencial. Me direis que para ese viaje no hacian falta alforjas, pero es bastante curioso porque eso nos lleva a

$$
(1 – \int_0^\infty d(t) e^{-rt} dt) = 1 – \int_0^\infty d(t) e^{-rt} dt = { r \over \lambda +r}
$$

con \(\lambda = 1/\tau_d\), y este seria el tiempo medio de deteccion, asi que $$R0 = k \tau_d$$.

La ecuacion de ayer para relacionar k con el exponente de crecimiento inicial se simplifica casi por sorpresa:
$$ k { r \over \lambda +r} = r e^{r t_l} $$
y ahora la solucion grafica es encajar dos una hiperbola \({ k \over \lambda +r}\) con la exponencial pura \(e^{r t_l}\) que condensa el tiempo de latencia. Este segundo termino vale siempre 1 en el origen, asi que la condicion de existencia de solucion es que la hiperbola debe ser mayor que uno en el origen, usease
$$ {k \over \lambda} > 1$$
¡que vuelve a ser la condición basica de la teoria epidemiologica: R0 > 1! Al menos esto parece que nos da bien.

Si no tuvieramos tiempo de incubacion la solucion de hecho seria bastante trivial:

$$r = k – {1 \over \tau_d}$$

Indicando que el exponente que vemos en las gráficas que salen por internet hay que esperar que tenga al menos dos contribuciones: la de la velocidad de contagio y la del tiempo de detección. Y por eso es por lo que tenemos siempre el debate sobre si queremos meterle presión a bajar una cosa o la otra.

Otro apunte. Si lo que queremos es precisamente ver la sensibilidad de r a cada uno de estos parametros puede que haya casos en los que en vez de resolver explicitamente (usando más aproximaciones) nos podamos apañar con los teoremas de derivadas parciales de funciones implicitas, usease consideramos

$$
0 = F(r(k,\lambda),k,\lambda) = {k \over \lambda +r} e^{rt_l}
$$

y aplicamos las reglas de encadenado aquellas de segundo de carrera

$$
{\partial r \over \partial k} = – {{\partial F \over \partial k} \over {\partial F \over \partial r} }
$$

etcetera… si conseguimos recordarlas y asegurarnos de que se cumplen las condiciones 😀