Programación de un modelo simple de epidemia

Los modelos son simplificaciones de la realidad (a veces verdaderas sobre-simplificaciones). Sería ilusorio pretender modelar toda la complejidad de un determinado problema; lo que sí podemos es modelar algunos aspectos, luego comprender -gracias al modelo- un poco más sobre el problema, y mejorar entonces el modelo, incorporando nuevos aspectos.

Un ente biológico tiene diferentes y numerosos aspectos; aquellos propiamente biológicos son bien complejos, no así los aspectos mecánicos. En estudios de genética, de crecimiento poblacional, de distribución espacial de organismos, de biología celular, etc. hay -en todos ellos- aspectos mecánicos, sencillos de simular y, a la vez, de gran interés práctico.
Aquí desarrollaremos un modelo de algunos aspectos mecánicos de una epidemia
(o de una pandemia, ya que ahora nos concierne tanto).

Consideremos una población, en un área determinada (una ciudad por ejemplo), con algunos individuos contagiados de un determinado virus (u otro ente transmisible); tal población, se verá afectada por nuevos contagios en virtud de su simple movilidad: tanto los individuos sanos como los infectados se mueven; de los contactos que ocurran, resultará nuevos contagios.

De lo dicho en este párrafo podemos enumerar cuatro parámetros importantes en la mecánica del contagio:

º El número de individuos o tamaño de la población
º El tamaño del área o hábitat de esa población
º El número inicial de contagiados dispersos en ese hábitat, y
º La relativa movilidad de los individuos. 

Esos parámetros serían entradas a nuestro modelo. La salida o resultados de interés podría ser una curva de número de infectados versus tiempo.

¿Cuáles han sido las mayores simplificaciones en este modelo?

No hemos considerado que algunos individuos pueden ser inmunes a la infección, que otros pueden padecer la infección por un tiempo y curarse (probablemente con ayuda médica), que algunos individuos no sobreviven la infección y por tanto decrecería el número de individuos en el tiempo considerado; por otra parte la población podría, por el contrario, crecer por nacimientos o por inmigración... Nuestro modelo dejará por fuera todas esas consideraciones.

¿Qué podemos hacer con un modelo tan simple? El modelo -concentrado en los aspectos mecánicos de la epidemia- puede ilustrar varios comportamientos de interés, como el efecto de la densidad de población (es decir de la relación entre el tamaño del hábitat y el tamaño de la población), del número inicial de infectados, de la movilidad, y quizá de otros parámetros fáciles de incorporar (como patrones de movilidad, clustering, aislamiento, confinamiento).

Desarrollo paso a paso un modelo general de epidemias

Modelar es asignar entes matemáticos o computacionales para representar entes de la realidad bajo estudio. El hábitat (realidad del problema) lo representaremos por una matriz (ente matemático). En cada elemento de la matriz puede haber o no un individuo, y a su vez, ese individuo puede estar sano o contagiado. Muchos elementos en esa matriz estarán vacíos (lo contrario sería un hacinamiento).

Necesitamos simular ahora que los individuos se mueven en esa matriz (en su hábitat). Consideraremos un movimiento aleatorio hacia sus vecindades (random walk). Claro que esto es otra simplificación, los desplazamientos reales no son random, son orientados, intencionales; sin embargo, si tal consideración es muy cierta a nivel de individuo, al mirar el conglomerado de toda una población lo que se observa es un ir y venir, cruzarse, devolverse, es decir algo bastante random. En ese deambular, sea intencionado o random, aparecen las oportunidades de contacto, y por tanto de contagio. Consideraremos (¡bendita sean las simplificaciones!) que cada contacto es una infección segura (asignar una probabilidad de infección distinto de uno no es nada complicado, de modo que no nos preocupemos con esta simplificación).

Quizá el lector acaba de descubrir cuál será la dinámica de esta simulación

En efecto, es ésa: partiremos de un número "pob" de individuos distribuidos aleatoriamente sobre un hábitat (matriz o grilla) "hab'. Dispersos también aleatoriamente, en ese hábitat, hay un número inicial de contagiados o enfermos "e1". A cada paso (simulation step) se recorre toda la matriz moviendo, cada individuo, hacia localidades vecinas; dicho movimiento puede ser nulo (no se mueve en este paso), o puede ser de uno, o dos, o más saltitos (según lo indique el parámetro movilidad, "mov") hacia arriba, abajo, derecha o izquierda. Aparecerán algunos problemas de borde (nunca faltan), pero no nos preocupemos de eso a estas alturas.

Si el movimiento aleatorio (el saltito) conduce a una celda vacía, simplemente se mueve a ella el individuo en consideración: podría haber estado sano o enfermo, igual se movió y está en una nueva posición en la grilla. Si el movimiento conduce en cambio hacia una posición ocupada, ocurre algo diferente: si el individuo en la posición de partida estaba sano, pero el individuo en la posición de llegada estaba enfermo, ha ocurrido un contagio; o viceversa, si el de la posición de partida estaba enfermo y el de llegada sano: igualmente ocurre un contagio (se incrementará el contador de contagiados o enfermos, "e"); otra posibilidad es que ambos individuos estaban ya enfermos, simplemente quedan, ambos, en esa misma condición (no se incremente "e"). Finalmente si ambos estaban sanos, siguen estando sanos. En cualquiera de estos casos, dado que no vamos a montar un individuo encima del otro, ambos quedan en sus posiciones, ya tendrán chance de moverse en los sucesivos pasos del algoritmo.

Cuando terminemos de recorrer toda la matriz o grilla hemos completado un paso y tenemos un nuevo total de contagiados, "e"; guardaremos ese valor actual en un vector que vamos a denominar "enfermos", para luego hacer una gráfica de número de enfermos versus tiempo.

En cada paso, el tiempo, "t", se incrementa en una unidad (la unidad que se decida en base a la naturaleza de la epidemia). Llamemos "lapso" el tiempo total (o tiempo final) de simulación. 

Consideraremos una grilla cuadrada de dimensión "dim", es decir nuestro hábitat es de tamaño dim x dim. Esta matriz podría ser de números enteros: un "0" representa elemento vacío, un "1" indica elemento ocupado, más exactamente indica que hay  un individuo sano allí; el valor "2" también indica que el elemento está ocupado, pero con un individuo enfermo. Nuestro programa resultará más legible si usamos nombres descriptivos para estos valores:

CONST libre = 0; sano = 1; enfermo = 2;

(Alternativamente la matriz "hab", en lugar de enteros, podría ser de records con dos flags: uno que indica ocupado / libre y el otro sano / enfermo. Nos hemos decidido por el 0, 1, 2).

Lo descrito se resume en el siguiente pseudocódigo

FOR t:= 0 TO lapso -1 DO
    Step (t);    --simulation step
    enfermos[t]:= e    --número de enfermos en el tiempo "t"
END;
MostrarCurva

Donde "Step" el algoritmo que implementa el paso es:

FOR i:= 0 TO dim -1 DO    --recorre filas y columnas, intenta mover c/elemento
   FOR j:= 0 TO dim -1 DO
      IF hay individuo en [i, j] THEN
         Next (i, j, x, y);    --escoger aleatoriamente posible localidad vecina
         IF hab[x, y] está ocupado THEN
            IF hab[i, j] # hab[x, y] THEN    --uno de ellos enfermo, el otro sano
               hab[i, j]:= enfermo; hab[x, y]:= enfermo;
               INC (e)
            END
            ELSE la localidad [x, y] está desocupada:
               hab[x, y]:= hab[i, j];    --mover [i, j] hacia esa posición
               hab[i, j]:= libre    --y marcar [i, j] libre
            END
        END
    END
END

Nos queda el detalle de escoger aleatoriamente la posible localidad vecina, [x, y], para intentar mover hacia ella el elemento  [i, j]. Quizá a primera vista no luce muy claro el siguiente trozo:

Next: --calcula [x, y] a partir de [i, j]
    x:= i + IO.Random (5) - 2;
    y:= j + IO.Random (5) - 2;
    IF (x <= 0) OR (x >= size) THEN x:= IO.Random (size) END;
    IF (y <= 0) OR (y >= size) THEN y:= IO.Random (size) END

Las dos primeras líneas invocan un generador de números pseudo-random. La función "IO.Random" genera, aleatoriamente un número entre 0 y el valor de su argumento disminuido en uno; en este caso, entre 0 y 4; como le restamos 2, estamos en definitiva generado, aleatoriamente alguno de estos números: -2, -1, 0, +1, +2. Ese pequeño desplazamiento se lo agregamos (o restamos según el signo) a las coordenadas actuales [i, j], para obtener la celda destino cercana [x, y]. Ahora debe estar más claro.

Las últimas dos líneas resuelven el problema de borde que habíamos anticipado. Aquí hemos decidido escoger algún lugar al azar en la grilla si la nueva posición quedase fuera de la grilla. Hay otras soluciones (como rotar, "por detrás de la grilla", enviado al extremo derecho si se trataba de salir por la izquierda, etc.).

Hagamos más general este trozo de código que llamamos "Next". Allí hay unas constantes (el 5 y el 2) que representan la movilidad de los individuos; para poder experimentar con ese parámetro, vamos a llamar "mov" al argumento que pasamos a Random, y que como vimos define el rango de números que devolverá esa función. Como el papel del -2 es mover ese rango para que ofrezca tanto valores negativos como positivos, lo cambiaremos por el cociente entero de dividir m / 2. La forma general sería entonces:

... IO.Random (mov) - mov DIV 2 --para lograr simetría respecto a cero, m debe ser impar

Ése es todo el programa de simulación. Sólo queda inicializar las variables, y disponer de funciones para hacer el gráfico y, eventualmente, para hacer un mapa de puntos con la distribución inicial y final de los individuos, o una animación con las sucesivas distribuciones que resultan en cada "Step".

Vamos a definir las variables (dim, pob, e1, etc.) y a colocar encabezamiento a cada uno de estos trozos de código para convertirlos en PROCEDURES. Todo ello constituye el módulo "Epid". 

Programa Final

Usaremos el lenguaje Component Pascal.

Una versión en Python se puede construir fácilmente a partir de este código.

MODULE Epid;   (* Simulador Simple de Epidemias *)
(*

Copyright (c) 2013 - 2017 BlackBox Framework Center
Copyright (c) 1994 - 2013 Oberon microsystems, Inc., Switzerland. All rights reserved.
Copyright (c) 2020 René Dorta Franceschi.

Disclaimer. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*)

IMPORT IO;

CONST
    libre = 0; sano = 1; enfermo = 2; maxTry = 10;

VAR
    dim*, pob*, e1*, lapso*, mov*, e, seed*: INTEGER;
    hab: POINTER TO ARRAY OF ARRAY OF INTEGER;
    legend: ARRAY 256 OF CHAR;
    
PROCEDURE Init;
    (* Llena el habitat con tantos individuos como diga "pob", entre ellos "e1" enfermos *)
    VAR n, i, j, try: INTEGER;
BEGIN
   IO.Assert (pob <= dim * dim, "¡Sobrepoblado!!");
   IO.Assert (e1 < pob, "Población ya infectada");
   IO.Assert (mov < dim, "Movilidad exagerada");    (* faltaría otras validaciones *)
   IF seed # 0 THEN IO.Seed (seed) END;
   IO.StrInt ("semilla RNG=", IO.seed); IO.Line;
   IF seed # 0
      THEN legend:= IO.C ("seed=",seed)$ ELSE legend:= IO.C ("seed=",IO.seed)$
   END;
   e:= 0;
   NEW (hab, dim, dim); (*  crea una grilla vacía (espacios libres) *)
   FOR n:= 1 TO pob DO
      try:= 0;
      REPEAT
         i:= IO.Random (dim); j:= IO.Random (dim);
         INC (try)    (* se reintenta por si el azar nos lleva a puesto ocupado *)
      UNTIL (hab[i, j] = libre) OR (try >= maxTry);
      hab[i, j]:= sano;
      IF e < e1 THEN
         hab[i, j]:= enfermo;
         INC (e);
      END
   END;
   legend:= legend + IO.C (" dim=",dim) + IO.C("x",dim) +
      IO.C (" pob=",pob) + IO.C (" mov=", mov) + IO.C (" e1=",e1);
   IO.StrLn (legend);
   IO.StrRealLn (" densidad=", pob / (dim * dim))
END Init;

PROCEDURE Next (i, j: INTEGER; OUT x, y: INTEGER);
BEGIN    (* devuelve [x, y], un vecino random de [i, j] *)
    x:= i + IO.Random (mov) - mov DIV 2;
    y:= j + IO.Random (mov) - mov DIV 2;
    IF (x <= 0) OR (x >= dim) THEN x:= IO.Random (dim) END;
    IF (y <= 0) OR (y >= dim) THEN y:= IO.Random (dim) END
END Next;

PROCEDURE Step (t: INTEGER);
   VAR i, j, x, y: INTEGER;
BEGIN
   FOR i:= 0 TO dim -1 DO    (* recorre filas y columnas, intentando mover c/elemento *)
      FOR j:= 0 TO dim -1 DO
         IF hab[i, j] # libre THEN    (* [i, j] es una posición ocupada *)
            Next (i, j, x, y);    (* escoger aleatoriamente posible localidad vecina *)
            IF hab[x, y] # libre THEN    (* [x, y] también es una posición ocupada *)
               IF hab[i, j] # hab[x, y] THEN    (* uno de ellos está enfermo, el otro sano *)
                  hab[i, j]:= enfermo; hab[x, y]:= enfermo;
                  INC (e)
               END
            ELSE  (* [x, y] desocupada: mover [i, j] hacia allí, y marcar [i, j] libre *)
               hab[x, y]:= hab[i, j];
               hab[i, j]:= libre
            END
         END
      END (*FOR j*)
   END (*FOR i*);
END Step;

PROCEDURE Sim*;    (* esto sería lo que otros lenguajes llaman "main" *)
   VAR t: INTEGER; enfermos: POINTER TO ARRAY OF REAL;
BEGIN
   Init;
   NEW (enfermos, lapso);
   FOR t:= 0 TO lapso -1 DO
      Step (t);
      enfermos[t]:= e;
      IF (t MOD (lapso DIV 10)) = 0 THEN IO.IntTab (t); IO.IntLn (e) END
   END;
   IO.Percent2 (e, pob); IO.Line;
   IO.Plot (enfermos, legend, "Infectados vs Tiempo", "", "")
END Sim;

BEGIN
   lapso:= 200; dim:= 100; pob:= 200; e1:= 5; mov:= 5  (* valores por defecto *)

END Epid.Sim

La interfaz con el usuario luce así (el botón "Sim" arranca la simulación):

Algunos resultados permiten ver cómo influye -en la rapidez de contagios- la densidad de población (cambiando dim y pob), el número inicial de infectados (cambiando e1, el numero inicial de individuos contagiados), la movilidad (cambiando mov) y los patrones de dispersión sobre el hábitat (cambiando la semilla del generador de números random).

Al lado uno de los típicos gráficos obtenidos --gracias a las rutinas para ploting en Component Pascal desarrolladas por Ivan Denisov. Con los parámetros mostrados sobre el propio gráfico, se observa la típica curva sigmoidea: muestra un crecimiento lento al comienzo (por el bajo número de infectados). A medida que ocurren nuevas infecciones el incremento de contagios se acelera, entrando a una sección lineal de la curva y luego se desacelera cuando se acerca a la saturación.

En la figura a arriba se compara el comportamiento de la infección en el tiempo para dos diferentes valores de movilidad (movilidad relativa 3 y 5). Se observa cómo ese factor incide notablemente en las pendientes de la curva. 

Para utilizar este software (experimentar con él, y extender el modelo) necesita copiar el módulo Epid (desde MODULE Epid hasta END Epid.),  y correrlo en el ambiente BlackBox que está disponible en el sitio web http://www.zinnamturm.eu/downloads.htm

¡A evitar contagios!

Vuestros comentarios y observaciones son bienvenidos.

Abril/2020