Introducción

Inspirado por este artículo de Washington Post: https://www.washingtonpost.com/graphics/2020/world/corona-simulator/ ; intenté replicar la simulación agregando algunos parámetros más, y la posibilidad de usar una forma distinta para la caja que se usa normalmente, para ilustrar la simulación con una forma conocida y crear más consciencia sobre cómo diferences deciones pueden cambiar los resultados de la situación actual.
Este es un modelo basado en agentes, basado en el caminante aleatorio para modelar el movimiento, y permite dar un polígono como “contenedor”.
Una diferencia importante con el modelo del artículo de Washington Post es que aquí estoy incluyendo también la posibilidad de muertos además de recuperados, y considerar además que existe una capacidad de atención limitada, lo que se refleja al considerar estados de gravedad con y sin cuidado.
En este notebook solo ejecuto una simulación para cada escenario, para obtener resultados estadísticos se deberían correr muchas simulaciones.

Especificaciones de la simulación

Estados

En este modelo se consideran estos estados:
◼
  • 0: Sano (suceptible de contagiarse)
    • ◼
  • 1: Infectado (incubando)
    • ◼
  • 2: Infectado con síntomas
    • ◼
  • 3: Infectado grave sin cuidados
    • ◼
  • 4: Infectado grave con cuidados
    • ◼
  • 10: Recuperado
    • ◼
  • 11: Muerto

    • Parámetros

    ◼
  • “ProbabilityContagion”: probabilidad de obtener la infección cuando se “toca” con un agente infectado.
    • ◼
  • “DaysAsymptomaticMin”: mínimo número de días para comenzar a mostrar síntomas.
    • ◼
  • “DaysAsymptomaticMax”: máximo número de días para comenzar a mostrar síntomas.
    • ◼
  • “ProbabilitySevere”: probabilidad de desarrollar un caso severo
    • ◼
  • “ProbabilitySevereDies”: probabilidad que un caso severo con cuidados muera
    • ◼
  • “ProbabilitySevereWithoutCareDies”: probabilidad que un caso severo sin cuidados muera
    • ◼
  • “HealthCapacity”: capacidad total para atender casos severos
    • ◼
  • “DaysOutcomeMin”: mínimo número de días para obtener un resultado (recuperado o muerto)
    • ◼
  • “DaysOutcomeMax”: máximo número de días para obtener un resultado (recuperado o muerto)
    • ◼
  • “TimeStep”: tiempo para cada paso usado para los movimientos
    • ◼
  • “ReinfectionPropability”: probabilidad que un agente recuperado se recaiga con la infección
    • ◼
  • “StatesInfecting”: alternativas con los estados que pueden infectar a otros. Por ejemplo, si todos los infectados infectan, o solo los que no tienen síntomas, o también los que no tienen cuidados y variantes así, que podrían simular aislamiento al 100% de casos detectados. Posibles valores recomendados: 1; 1 | 2; 1 | 2 | 3; 1 | 2 | 3 | 4.
    • Los parámetros Min/Max son usados para calcular una distribución normal alrededor de su media y con desviación estándar de
    1
    max-min
    .

    Código de la simulación

    Estructuras de datos

    Usaremos una lista de posiciones enteras, y una lista de estados de la forma {estado, días infectado, movimiento}, con los siguientes valores:
    ◼
  • estado: uno de los valores de estado {0, 1, 2, 3, 4, 10, 11}
    • ◼
  • días infectado: contador de tiempo usado solo en los casos infectados
    • ◼
  • movimiento: si un agente se mueve es 1; si no se mueve es 0; y es -1 si se movía pero ahora es detectado como infectado y ya no se mueve (pero se moverá cuando se recupere)
    • Finalmente, la estructura que describirá el sistema será dada por {posiciones, estados}.

    Constantes

    Asociemos un color con cada estado:
    In[]:=
    statecolor=Association[0Gray,1Pink,2Red,3Darker@Red,4Brown,10Green,11Black]
    Out[]=
    0
    ,1
    ,2
    ,3
    ,4
    ,10
    ,11
    
    Patrón con los estados infectados:
    In[]:=
    $infectedstates=1|2|3|4;

    Creación de la simulación

    Queremos generar las condiciones iniciales para un contenedor dado (definido por una región, entidad, polígono). Reescalaremos la región pues queremos que las el número de agentes tenga cierta densidad objetivo.
    Data la región o entidad y el número de agentes, retornaremos la región reescalada, una función de pertenencia a la región, y las posiciones iniciales con sus estados.

    Sobre la región

    Necesitamos una función que modifique la región a utilizar y la reescale de tal manera que los agentes puedan entrar usando sus posiciones expresadas en números enteros, y así asegurarnos que haya espacio para ellos. Esta función retornará la región modificada, la función de pertenencia a la región correspondiente y la función de transformación de escala usada para pasar de la región original a la nueva. Esa función de transformación será usada luego para transformar las coordenadas usadas como referencia.
    In[]:=
    Options[prepareRegion]={};​​prepareRegion[r:_Region|_Entity|_Polygon,targetarea_]:=Module[{region,area,newregion},​​If[MatchQ[r,_Entity],​​region=EntityValue[r,"Polygon"]/.g_GeoPositiong["LongitudeLatitude"];​​,​​region=r;​​];​​area=Area[region];​​If[!NumericQ[area],​​Return[{$Failed,$Failed,$Failed}];​​];​​newregion=RegionResize[region,Abs[Subtract@@RegionBounds[region][[1]]]Sqrt[targetarea/area]];​​{newregion,RegionMember[newregion],RescalingTransform[RegionBounds[region],RegionBounds[newregion]]}​​];

    Condiciones iniciales

    Las condiciones iniciales para la simulación son las posiciones de los agentes y sus estados. Para generar los puntos aleatorios dentro de la región usaremos RandomPoint, pero al redondearlos luego, necesitamos usar Select para verificar que los agentes estén dentro y borrar duplicados; después de eso necesitamos agregar más de forma iterativa, para asegurarnos obtener la cantidad requerida (y siempre agregando un contador en el bucle para evitar que se ejecute infinitamente).
    In[]:=
    createPositions[n_,region_,regionq_:Automatic]:=Module[{memberq=Replace[regionq,AutomaticRegionMember[region]],positions,i=0},​​positions={};​​While[Length[positions]<n&&i++<n,​​positions=DeleteDuplicates[Join[positions,Select[Round@RandomPoint[region,n-Length[positions]],memberq]]];​​];​​If[Length[positions]<n,​​$Failed,​​positions​​]​​];
    Para inicializar el estado de cada agente usamos un valor aleatorio de 1 o 0, en el que se usa la proporción de aislamiento para definir si se mueven o no. Luego se usarán distintas estrategias para establecer los infectados. Puede ser un número especificado de infectados (se escogerán cuáles son al azar), o podemos dar una lista establecida (en caso se calcule en otra función externa), o podemos usar UpTo como cabecera para especificar que queremos los primeros índices como infectados.
    In[]:=
    createStates[n_,infected:{__Integer}|_Integer|UpTo[_Integer],isolationProp_]:=Module[{states},​​states=Table[{0,0,RandomChoice[{isolationProp,1-isolationProp}{0,1}]},n];​​If[!TrueQ[And@@MatchQ[infected,_List]],Return[$Failed]];​​Switch[infected,​​_List,​​states[[infected]]=ConstantArray[{1,0,1},Length[infected]];,​​_UpTo,​​states[[;;infected[[1]]]]=ConstantArray[{1,0,1},infected[[1]]];,​​_,​​states[[RandomSample[Range[n],infected]]]=ConstantArray[{1,0,1},infected];​​];​​states​​];
    In[]:=
    searchInfectionFocus[posisitions_,point_,transformation_,i_]:=Nearest[posisitions"Index",transformation[point],i];

    Simulación

    La función principal para generar el estado inicial de la simulación llama a las funciones de inicialización previas. Esta función es la que retorna la región, la función de pertenencia, las posiciones iniciales y estados.
    In[]:=
    Options[createSimulation]:={"Infected"1,"Density"Automatic,"InfectionPoint"None,"Isolation"0.1}​​createSimulation[n_,r:_Region|_Entity|_Polygon,OptionsPattern[]]:=Module[{region,regionq,positions,transform,infected,density=OptionValue["Density"]/.Automatic1/10,ninfected,isolation,states},​​{region,regionq,transform}=prepareRegion[r,n/density];​​positions=createPositions[n,region];​​ninfected=OptionValue["Infected"];​​Which[MatchQ[OptionValue["InfectionPoint"],_Entity],​​infected=searchInfectionFocus[positions,EntityValue[OptionValue["InfectionPoint"],"Position"]["LongitudeLatitude"],transform,ninfected];​​,​​MatchQ[OptionValue["InfectionPoint"],{_?NumericQ,_?NumericQ}],​​infected=searchInfectionFocus[positions,OptionValue["InfectionPoint"],transform,ninfected];​​,​​True,​​infected=ninfected;​​];​​isolation=OptionValue["Isolation"];​​states=createStates[n,infected,isolation];​​{region,regionq,positions,states}​​];

    Evolución

    Rutinas principales

    La rutinas principales son las encargadas ejecutar la simulación completa y calcular la evolución en cada paso.
    Para cada paso hacemos 3 cosas:
    ◼
  • Primero: calculamos el movimiento de los agentes, evitando sobreponerlos en la misma posición y verificando que no se salgan de la región especificada (para lo que se usa la función de pertenencia a la región).
    • ◼
  • Segundo: calculamos las nuevas infecciones e incrementamos el contador de tiempo de la infección de los agentes ya infectados. En esta etapa es que se calcula si un agente se pone peor, se recupera o muere.
    • ◼
  • Tercero: calculamos cuáles obtienen atención, lo cual se basa en la capacidad total de atención. Esto se hace considerando la capacidad total, y cuántos están en el estado “infectados con cuidados”, y si hubiera disponibilidad, algunos casos severos sin cuidados se cambiarán a casos con cuidados.

    • Funciones auxiliares

    Las funciones auxliares son las que se llaman desde las rutinas principales.
    Para calcular cuál agente es vecino de cuál, primeros los ordenamos por su coordenada x, luego por la coordenada y, y después de agruparlos según estas coordenadas, tenemos un grupo más pequeño en el que podemos verificar cuáles son realmente vecinos.
    Esta es una de las rutinas centrales, cómo un agente infectado evoluciona. Se toma en consideración todos los parámetros de probabilidades y tiempos. De acuerdo a valores aleatorios, un agente infectado puede mostrar los síntomas en algún caso, pero según los pesos de probabilidades puede ser un caso leve o severo, o si tuviera mucha suerte, incluso recuperarse directamente. Si es un caso ligero, eventualmente se recuperará, y solo en los casos severos es que podrán haber recuperación o muerte (dependiendo esa probabilidad según sea un caso con cuidados o no).

    Otras funciones auxiliares

    Estas funciones auxiliares ayudan a determinar si se ha alcanzado un punto estable en la ejecución de una simulación, esto es, que no haya infectados.

    Funciones de visualización

    Creamos una función principal de visualización, la cuál recibe todos los pasos de la ejecución de la simulación, y la región, y puede graficarlo para cualquier paso intermedio.

    Funciones auxiliares de visualización

    Ejecutando un caso de simulación

    Simulación con Perú como fondo

    Necesitamos definir los parámetros para la simulación en una asociación (Association). Establecemos una probabilidad de mortalidad de 10% de los casos severos con cuidados (que ya son un 20%, lo que da un global de 1%), y de 80% para los casos severos sin cuidados. También establecemos el paso de tiempo como 1h. Otro parámetro importante es la capacidad de cuidados, que estableceremos como 3 para esta población total de 200.

    Simular con 10% de aislamiento

    Para inicializar la simulación, usamos la entidad Perú, y establecemos además que el primer infectado se localice cerca de Lima. Debemos recordar que este modelo está usando el polígono como contenedor, y no debe confundirse con un modelo predictivo que esté usando un área geográfica.
    Ejecutamos la simulación por 90 días.
    Para crear una visualización interactiva más ligera (con Manipulate) para cada paso, usaremos como paso final el paso en que ya no haya infectados.

    Ahora con 90% de aislamiento

    Para simular nuevamente con solo el parámetro de aislamiento distinto, no necesitamos retornar la región ni función de pertenencia, pues ya los tenemos del caso anterior. Entonces ahora ejecutamos y creamos la visualización para este caso.

    Visualizar ambos casos juntos

    Podemos poner ahora ambas simulaciones, que se diferencian solo por un parámetro, en la misma visualización interactiva, y se usará como paso final mayor tiempo de estabilidad.

    Simulación con la isla de Gran Bretaña de fondo

    Para Gran Bretaña usaremos los mismos parámetros, pero incrementaremos la capacidad de atención.

    Polígono de la isla de Gran Bretaña

    Si queremos usar una proyección bonita, no podemos usar la extracción automática usada en el código, pero podemos proyectar el polígono manualmente:

    Simular con 10% de aislamiento

    Calculamos ahora las condiciones iniciales para este caso, y ejecutamos la simulación por 90 días.

    Ahora con 90% de aislamiento

    Cuando cambiamos al 90%, no necesitamos retornar la región ni función pertenencia.

    Visualizar ambos casos juntos

    Ahora pondremos ambas simulaciones juntas, como se hizo para el caso de Perú antes.

    Notas finales

    Este es un modelo sencillo donde no se intenta imitar la realidad, pero intenta ilustrar cómo distintos comportamientos pueden cambiar el resultado final.
    El hecho que se use un polígono como contenedor no cambia demasiado el resultado final, pero ayuda a visualizar el impacto en casa.
    Podríamos tratar de incorporar más parámetros, como agregar otro parámetro para la probabilidad de contagio de los agentes estáticos. Otro parámetro interesante podría ser tener un aislamiento variable en el tiempo, que podría utilizar para simular casos en los que la decisión de aislamiento es tomada tarde, o si se quitara después de cierto período, pero antes que la infección hubiera desaparecido.
    Una versión preliminar de este código es la que se usó para crear esta animación: https://twitter.com/fjsistemas/status/1239390168848707585
    ¡Quédense en casa!