No soy epidemiólogo, pero imparto clases de métodos matemáticos para ingenieros. En autocuarentena en casa, he decidido implementar en Matlab un modelo sencillo de la epidemia de coronavirus. En lugar de los conocidos modelos SIR y SEIR, me he decantado por el modelo SEIRV con tasas de transmisión variables, que es un poquito más realista. Por supuesto, estimar el valor de sus parámetros es difícil, así que he tomado los valores que se acaban de publicar para la epidemia en Wuhan, China. Quizás sea de interés para algunos de los lectores de este blog. Pero, cuidado, siendo un modelo matemático tan sencillo, que tiene en cuenta tan pocos factores, sus resultados no deben ser extrapolados sin ton ni son. La epidemiología rigurosa requiere usar modelos más complejos que incluyan toda la información disponible.
El modelo de juguete que voy a presentar es solo eso, un juguete que me servirá para ilustrar el uso de Matlab a la hora de resolver problemas de valores iniciales en ecuaciones diferenciales ordinarias. Por ello me centraré en lo fácil que es programar en este lenguaje este tipo de modelos matemáticos. He seleccionado el lenguaje Matlab porque es el que uso en mis clases para ingenieros (uso R en mis clases para bioquímicos); quienes prefieran una alternativa
open source pueden usar Octave de GNU, que es un lenguaje casi idéntico. El código se puede adaptar de forma sencilla a otros lenguajes matemáticos, por ejemplo, a Mathematica usando NDSolve; pero en otros lenguajes hay que usar librerías o paquetes específicos, por ejemplo, deSolve en R, o diffeqpy en Python.
El artículo (
open access) que me ha servido en bandeja el modelo y sus parámetros es Chayu Yang, Jin Wang, “A mathematical model for the novel coronavirus epidemic in Wuhan, China,” Mathematical Biosciences and Engineering 17: 2708-2724 (11 Mar 2020), doi:
https://doi.org/10.3934/mbe.2020148.
% componentes del vector y = SEIRV
S=1;E=2;I=3;R=4;V=5;
dSEIRVdt = @(t,y) [
% dS/dt =
Lambda-betaE(y(E))*y(S)*y(E)-betaI(y(I))*y(S)*y(I)-betaV(y(V))*y(S)*y(V)-mu*y(S)
% dE/dt =
betaE(y(E))*y(S)*y(E)+betaI(y(I))*y(S)*y(I)+betaV(y(V))*y(S)*y(V)-(alfa+mu)*y(E)
% dI/dt =
alfa*y(E)-(w+gamma+mu)*y(I)
% dR/dt =
gamma*y(I)-mu*y(R)
% dV/dt =
xi1*y(E)+xi2*y(I)-sigma*y(V)
];
Introducir las ecuaciones diferenciales dY/dt = F(t,Y) en Matlab es muy fácil, basta definir la función F(t,Y). En este caso tenemos cinco ecuaciones acopladas, luego Y(t) será un vector de cinco componentes. He llamado a la función dSEIRVdt(t,Y) donde Y(1) es S, Y(2) es E, etc. En este código he usado una función anónima @(t,y) en Matlab. A mis estudiantes les recomiendo que no abusen de las funciones anónimas, porque el resultado es un código mucho más difícil de entender que cuando se usan funciones. Te muestro la diferencia.
function dYdt = dSEIRVdt (t,SEIRV)
% function dYdt = dSEIRVdt (t,SEIRV) % modelo SEIRV que usa constantes globales
%%
global Lambda betaE betaI betaV mu alfa w gamma xi1 xi2 sigma;
S=SEIRV(1); E=SEIRV(2); I=SEIRV(3); R=SEIRV(4); V=SEIRV(5);
%%
dSdt = Lambda-betaE(E)*S*E-betaI(I)*S*I-betaV(V)*S*V-mu*S;
dEdt = betaE(E)*S*E+betaI(I)*S*I+betaV(V)*S*V-(alfa+mu)*E;
dIdt = alfa*E-(w+gamma+mu)*I;
dRdt = gamma*I-mu*R;
dVdt = xi1*E+xi2*I-sigma*V;
%%
dYdt = [dSdt ; dEdt; dIdt; dRdt; dVdt ]
end
En el modelo SEIRV se usan cuatro poblaciones de humanos: susceptibles (S) son las personas que pueden ser infectadas; expuestas (E) al virus son las personas que lo contagian, pero que no presentan síntomas durante el periodo de incubación, por lo que no requieren hospitalización y aislamiento riguroso; infectadas (I) son las personas que presentan síntomas y están hospitalizadas, sin contacto con las expuestas; y recuperadas (R) son las personas que se han inmunizado contra el virus tras superar la infección; además, se incluye la concentración de virus (V) en el entorno debida a las personas expuestas (pues las infectadas influyen poco al estar hospitalizadas, sin contacto directo con las susceptibles).
Un modelo tan sencillo tiene 12 parámetros (los 11 mostradosmás un parámetro adicional común a todas las funciones beta) cuyos valores deben ser ajustados con datos epidemiológicos. Un ajuste completo es inviable (explorar un espacio con 12 dimensiones con 5 muestras por dimensión requiere unos 250 mil puntos; pero aún no hay tantos infectados). En la práctica, se estiman algunos parámetros de forma razonable y razonada, y se ajustan unos poquitos parámetros.
El artículo en que me baso estima 9 parámetros a partir de la literatura previa y ajusta 3 parámetros usando los datos epidemiológicos de Wuhan entre el 23 de enero de 2020 y el 10 de febrero de 2020, es decir, en la fase de crecimiento exponencial. Por tanto, dicho ajuste debe ser tomado como un simple ejemplo; en ningún caso será realista si se extrapola a otros países o a toda la pandemia.
% parametros del modelo SEIRV
http://doi.org/10.3934/mbe.2020148
Lambda = 271.23 ;% flujo de poblacion hacia Wuhan (por dia)
betaE0 = 3.11e-8;% tasa de transmision entre S y E (/personas/dia)
betaI0 = 0.62e-8;% tasa de transmision entre S e I (/personas/dia)
betaV0 = 1.03e-8;% tasa de transmision entre S y V (/personas/dia)
c = 1.01e-4;% coeficiente de ajuste de la transmision
mu = 3.01e-5;% tasa natural de fallecimientos (por día)
alfa = 1/7 ;% 1/alfa = periodo de incubacion (/dias)
w = 0.01 ;% tasa de fallceimientos por la infeccion (por dia)
gamma = 1/15 ;% tasa de recuperacion de la infeccion (por dia)
sigma = 1 ;% tasa de eliminacion del virus del ambiente (por dia)
xi1 = 2.30 ;% dispersion del virus por infectados asintomaticos
xi2 = 0 ;% y por infectados aislados (por persona por dia por ml)
El parámetro Λ representa el flujo de entrada de personas a Wuhan; μ es la mortalidad de la población sin relación con la infeccion; 1/α es el periodo de incubación en el que la infección es asintomática (7 días); w es la mortalidad debida a la infección; 1/γ es el número de días para la recuperación de un infectado (15 días); ξ1 y ξ2 es la contribución a la concentración ambiental del virus debida a las personas expuestas e infectadas, resp.; y σ es el ritmo de degradación del virus en el entorno (tanto por causas ambientales, como por desinfecciones programadas). Las funciones βE(E) y βI(I) representan las tasas de contagio de persona a persona entre las expuestas y las susceptibles, y entre las infectadas y las susceptibles, resp.; la función βV(V) representa la tasa de contagio de una persona debida a la concentración de virus en el entorno. Estas tres funciones decrecen de forma monótona, pues a mayores valores de E, I y V se tomarán medidas de control más drásticas para reducir el contagio.
% funciones de contagio (positivas y decrecientes con su argumento)
betaE = @(E) betaE0/(1+c*E); % para expuestos
betaI = @(I) betaI0/(1+c*I); % para infectados
betaV = @(V) betaV0/(1+c*V); % para el ambiente
En el artículo se usa la misma expresión matemática para las tres funciones que describen el contagio; para no introducir tres nuevos parámetros, solo se usa un único parámetro c común a las tres, que depende de las medidas de control que tratan de evitar los contagios. Para c=0 se obtiene el modelo SEIRV con parámetros constantes.
Esta figura del artículo muestra que para c=0 los resultados no describen la evolución de la epidemia en Wuhan desde el 23 de enero de 2020 (día 0 en la figura). El modelo SEIRV con parámetros constantes predice casi tres millones de contagiados y todos sabemos que en Wuhan ni siquiera se han alcanzado los cien mil. Esto nos recuerda que ajustar el modelo con parámetros estimados en fecha muy temprana (antes del 10 de febrero) conlleva una predicción muy exagerada del pico máximo de infectados. En redes sociales mucha gente usa una simple función exponencial con el número de reproducción y obtienen valores anuméricos sin ningún sentido. Ten mucho cuidado cuando alguien te muestre ese tipo de cálculos; son simplemente basura anumérica.
Esta figura muestra la predicción del modelo para c estimado usando los datos entre 23 de enero y 10 de febrero. Predice un pico de personas infectadas de unas 50 mil, lo que es mucho más realista. Sin embargo, predice una caída del número de personas infectadas mucho más lenta de la observada en la epidemia en Wuhan (300 días es casi un año). Además, predice un número de personas expuestas (infectadas asintomáticas) inferior a las que acaban siendo infectadas. Habrás leído en redes sociales a mucha gente estimar que tiene que haber entre 10 y 30 veces más personas expuestas que infectadas; este modelo no da cuenta de ese hecho (incluso aunque inicialmente haya más personas expuestas que infectadas).
El modelo SEIRV del artículo es muy sencillo. Sus conclusiones deben ser tomadas con mucha cautela. El artículo me ha gustado porque incluye un cálculo de los puntos de equilibrio y demuestra que para una número de reproducción R0>1 siempre existe un punto de equilibrio no trivial. Esta figura te muestra el plano de fases para las personas expuestas (E) e infectadas (I). El punto fijo calculado para R0 = 4.25 es E* = 1353, e I* = 2735, con S* = 2 583 683, R* = 6 528 015, y V* = 3111. Este punto fijo implica que la infección se vuelve endémica; siempre habrá un reservorio de virus entre la población. Aunque el objetivo de las medidas de la OMS y de muchos Gobiernos ha sido impedir este hecho, cada día parece más difícil evitar que así ocurra con la pandemia de COVID-19. Pero, repito de nuevo, el modelo SEIRV es muy sencillo y sus conclusiones debe ser puestas en cuarentena. Una vez pase la pandemia, gracias a modelos mucho más complicados, se podrá saber si el nuevo coronavirus se ha convertido en un nuevo coronavirus humano estacional.
% condicion inicial
S0 = 8998505 ;% susceptibles de infectar
E0 = 1000 ;% expuestos (infectados asintomaticos)
I0 = 475 ;% infectados (con sintomas y aislados)
R0 = 10 ;% recuperados
V0 = 10000 ;% concentracion de virus en el ambiente
Volvamos a mi código en Matlab. La condición inicial que se usa en el artículo está basada en información oficial de Wuhan (
fuente). Se considera una población de susceptibles de casi 9 millones de habitantes en Wuhan, 475 personas infectadas, 10 recuperadas y 1000 personas expuestas, todo ello a fecha de 23 de enero de 2020. La carga viral inicial en el entorno es muy difícil de estimar, pero los autores proponen un valor de 10 mil partículas virales por mililitro; no hay ningún razonamiento en el artículo que justifique este número, que me parece que a sido elegido más a ojo que otra cosa.
% resolucion numerica
SEIRV0 = [S0; E0; I0; R0; V0];
[t,SEIRV] = ode23t(dSEIRVdt,[0,3000],SEIRV0);
La ventaja de usar Matlab (o lenguajes matemáticos como Mathematica, R, Python, etc.) es que la resolución del sistema de ecuaciones diferenciales ordinarias solo requiere una única línea. En este código uso la función ode23t, que se suele usar en problemas moderadamente rígidos (moderate stiffness). ¿Por qué la uso? Porque he probado con todas…
% probar toda la odesuite
odesuite = {“ode23″,”ode45″,”ode15s”,”ode23s”,”ode23t”,”ode23tb”};
for ii=1:length(odesuite)
time = cputime;
eval(strcat(‘[t,SEIRV]=’,odesuite{ii},'(dSEIRVdt,[0,30000],SEIRV0);’));
cost(ii) = cputime-time,
end
He probado con todos los métodos de la odesuite de Matlab. El coste en tiempo de CPU ha sido de 0.58, 1.0, 0.047, 0.016, 0.016, y 0.016 segundos. Así que los métodos ode23s, ode23t, y ode23tb tardan unas 66 veces menos que ode45; por ello he elegido ode23t.
clf;
SEIRV0 = [S0; E0; I0; R0; V0];
[t,SEIRV] = ode23t(dSEIRVdt,[0,300],SEIRV0);
subplot(2,2,1); plot(t,SEIRV
,E),t,SEIRV
,I)); title(‘E,I’); xlabel(‘dias’);
[t,SEIRV] = ode23t(dSEIRVdt,[0,3000],SEIRV0);
subplot(2,2,2); plot(t,SEIRV
,E),t,SEIRV
,I)); title(‘E,I’); xlabel(‘dias’);
Este código y esta figura muestra los resultados de Matlab para el modelo SEIRV. El pico de la curva es razonable pero la caída del número de personas infectadas ha sido mucho más rápida (300 días es casi un año). Sin lugar a dudas, el método requiere una nueva estimación de sus parámetros.
% dibuja el plano de fase
I0v = [ 5 10 15 15 15 15 15 10 6 3 0 0 0 0 ]*1e4;
E0v = [ 0 0 0 1 2 3 4 4 4 4 3.5 2.5 1.5 0.5 ]*1e4;
clf; box on;
for iv=1:length(I0v)
SEIRV0 = [S0; E0v(iv); I0v(iv); R0; V0];
[t,SEIRV] = ode23t(dSEIRVdt,[0,50000],SEIRV0);
subplot(2,2,1); plot(SEIRV
,E),SEIRV
,I),’-‘); hold on;
subplot(2,2,2); semilogx(t,SEIRV
,E),t,SEIRV
,I)); hold on;
end;
subplot(2,2,1); xlabel(‘E’); ylabel(‘I’); title(‘phase portrait’);
subplot(2,2,2); xlabel(‘t’); title(‘E, I’); axis([0.1 5e4 0 1.5e5]);
Este código y esta figura muestran el plano de fases (E, I) a la derecha, así como las soluciones E(t) e I(t) en función del tiempo (en días). Se observa que las soluciones tienden a un punto intermedio alrededor de t = 100 días con E(100) ≈ 21 000, e I(100) = 39 000; tras este punto, todas las soluciones en el pano de fase evolucionan hacia el mismo punto de equilibrio E* = 1347, e I* = 2509. Este punto no coincide con el que aparece en el artículo, pero se puede calcular directamente en Matlab.
% calcula punto de equilibrio
SEIRVstar = fsolve(@(y) dSEIRVdt(0,y), SEIRV0)
SEIRVstar([E, I])
Recuerda que un punto de equilibrio del sistema de ecuaciones diferenciales ordinarias dY/dt = F(Y) es un punto Y* tal que F(Y*)=0; por tanto dicho punto es estacionario (usado como condición inicial el sistema se mantiene en ese punto como solución constante). El cálculo de Matlab ofrece el mismo resultado que se ha determinado con ode23t, que no coincide con el artículo. Supongo que la razón es el error introducido por el método numérico usado por los autores (que no aclaran cuál es en su artículo).
En resumen, sirva esta pieza como ejemplo del uso de Matlab para resolver ecuaciones diferenciales ordinarias. Los resultados epidemiológicos mostrados son poco rigurosos, pues el modelo SEIRV es muy sencillo (y yo creo que habría que reestimar sus parámetros a partir de los resultados actuales de la epidemia en Wuhan para obtener un mejor ajuste). Sin embargo, es un modelo mucho mejor que usar una ley exponencial, aunque la hayamos visto tantas veces en redes sociales que duele a los ojos tanto anumerismo. Si conoces Matlab (u Octave) , te recomiendo ejecutar el código con diferentes parámetros; usa los comentarios si logras resultados que merezca la pena remarcar.
y no nos pudó mostrar más de la realidad de China.
