Vamos a ver un ejemplo de lo fácil que es realizar algunas simulaciones simples usando autómatas celulares. En este caso modelaremos la evolución de una enfermedad en una población. Hay modelos matemáticos muy elaborados sobre enfermedades y que desde luego son más realistas que el que vamos a construir (aunque realmente se puede elaborar tanto como quieras).  La ventaja de los autómatas celulares es que son mucho más visuales lo que permite entender el desarrollo del sistema de forma más intuitiva, además sus mecanismos son más fáciles de entender y modificar que un sistema de ecuaciones diferenciales.

Para modelar nuestro autómata celular vamos a usar terra.js

Usaremos un autómata con 8 vecinos, cada celda va a simular ser un individuo diferente en contacto con otros individuos (sus celdas vecinas). Un individuo puede estar en cuatro estados:

• Susceptible: el individuo no ha sido infectado, pero puede contagiarse si cualquiera de sus vecinos está enfermo
• Infectado: el individuo esta enfermo, durante unos turnos seguirá en ese estado y será contagioso. Pasados esos turnos pasará a estar recuperado
• Recuperado: tras pasar la infección el individuo ni es contagioso ni puede volverse a contagiar
• Inmune: el individuo es inmune a la infección, no puede estar enfermo. Para el caso es como si fuera un recuperado solo que desde antes de infectarse, Sirve para simular inmunes y vacunados.

Esto es lo que se conoce como «Modelo SIR» (Susceptible Infectado Recuperado).

Parámetros

Tendremos que configurar varios parámetros para modelar la enfermedad.

Lo días (ciclos) que dura la enfermedad, illDuration, durante esos días el individuo infectado puede infectar a sus vecinos

Para caracterizar lo «infeccioso» que es el virus usaremos el número reproductivo básico o R0. Podemos interpretar este número como cuantas infecciones causa cada infectado en el periodo que está enfermo. En nuestro autómata celular cada ciclo es un día por lo que el número de individuos infectados por día será (en teoría) R0/illDuration. Necesitamos traducir esto a probabilidades que tiene cada uno de los 8 vecinos de contagiarse Pc = (R0/illDuration)/8. ¿Y qué pasa si R0 es mayor que 8? Es cierto que ningún individuo podrá infectar a más de 8 vecinos, este R= mayor se traducirá en una propagación más rápida de la enfermedad.

config.illDuration = 14; //duration of illness (cicles)
config.pc = (config.r0/config.illDuration)/8;

Cuando un individuo sano tiene como vecino a uno infectado tiene una probabilidad Pc de contagiarse. Esta probabilidad aumenta con el número de vecinos. La probabilidad de infectarse es 1-(1-Pc)^n siendo n el numero de vecinos infectados. ¿De donde sale esa formula?. Vamos despacito, si tenemos solo un vecino la probabilidad de infectarse es: Pc. Pero si tenemos dos vecinos la probabilidad de infectarse es: la probabilidad de infectarse del vecino A, la probabilidad de infectarse del vecino B y la posibilidad de infectarse de ambos. Resulta más sencillo calcular lo contrario, la probabilidad de no-contagiarse de ninguno, que para un vecino es Pnc = 1-Pc. La probabilidad de no infectarse de 2 vecinos es Pnc = (1-Pc)(1-Pc). De tres: Pnc = (1-Pc)*(1-Pc)*(1-Pc). De n: Pnc = (1-Pc)^n. Como lo que queremos saber es la probabilidad de infectarse y lo que hemos calculado es la de no infectarse, la de infectarse debe de ser: 1- Pnc = 1-(1-Pc)^n . Cada turno, cada individuo sano se genera un número al azar entre 0 y 1, si es menor que la probabilidad de infectarse de sus vecinos su estado cambia a infectado.

process: function (neighbors, x, y) {    
  if(!this.isIll && !this.wasIll && !this.isInmune){ //individuo sano 
    //vecinos enfermos alrededor
    var surrounding = neighborsCount(neighbors, 'isIll', true);
    if(Math.random() < 1-Math.pow(1-config.pc, surrounding)){
      this.isIll = true;          
      this.timesIll = config.illDuration;
    }
  } else if(this.isIll){ //individuo enfermo
      ......
  }
  return true;
}

Si un individuo ya esta enfermo le restamos tiempo de enfermedad cada turno hasta que llegue a 0 que es cuando deja de estar enfermo, deja de ser contagioso y pasa a estar recuperado.

process: function (neighbors, x, y) {    
  if(!this.isIll && !this.wasIll && !this.isInmune){ //individuo sano 
    ......
  } else if(this.isIll){ //individuo enfermo
    this.timesIll--;
    if(this.timesIll == 0){
      this.isIll = false;
      this.wasIll = true;
    }
  }
  return true;
}

Para caracterizar la población se puede establecer el tamaño, que porcentaje comienza enfermo y que porcentaje es inmune/vacunado. Al principio de la simulación los individuos se generan aleatoriamente según las probabilidades indicadas:

this.isIll = withProbability(config.startIll/100, true, false);
if(this.isIll){ //enfermo
  this.timesIll = config.illDuration;
} else { 
  this.timesIll = 0;
   //inmune
  this.isInmune = withProbability(config.startInmune/100, true, false);
}
this.wasIll = false;

Indicadores

A parte del tablero del autómata celular contamos con varios indicadores:

El numero de enfermos en ese ciclo

El número de recuperados hasta ese momento

Rt, el valor R0 del que hemos hablado antes solo es valido en el ciclo 0 (de hay el 0 detrás de la R) según la población se va infectando y recuperando el número decrece ya que hay menos población a la que infectar. Rt aproxima el valor de R en ese ciclo se calcula como: Rt = R0 * (1 – (enfermos+recuperados)/población). Cuando Rt es menor que 1 indica que la propagación empieza a remitir, con Rt mayor que 1 la propagación de la enfermedad aumenta.

Inmunidad de grupo (herdImmunity), es el punto en que Rt será menor que 1. Como ya hemos visto a partir de la cantidad de infectados y recuperados se puede estimar que valor de Rt tendrás. Por lo tanto podemos decir Rt = 1 = R0 * herdImmunity => herdImmunity = (R0-1)/R0;

El porcentaje de gente que nunca ha sido infectada, indica el avance de la enfermedad y su extensión.

stats.ills = ills;
stats.noInfectedPer = 100 * (config.population - (ills + recovered))/config.population;
stats.recovered = recovered;
stats.rt = config.r0 * (1 - ((ills+recovered)/config.population))
stats.herdImmunity = (config.r0-1)/config.r0;

El código completo puede verse aquí y una demostración aquí

Puedes ver un vídeo explicativo en mi canal de Youtube: