Modelización de la evolución de epidemias (grupo A14)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Modelo para epidemias. Grupo A-14 |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2014-15 |
| Autores | Pablo Busca, Manuel Escudero, Ángel Díaz, Miguel Ángel Cubo, Pablo Bueno |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
En este artículo se recogen algunos modelos sencillos de la evolución a lo largo del tiempo de una enfermedad infecciosa. Se supone que la epidemia se expande por una ciudad de N habitantes.
1 Modelo logístico
El crecimiento de una epidemia puede asimilarse al crecimiento de una población mediante el modelo logístico en cuanto a que se mide el crecimiento de un número de infectados con el límite que impone que dicho número de infectados no puede superar el de la población total N.
Así, se supone que la enfermedad es incurable y que todos los individuos infectados son contagiosos. Se llama I(t) al número de infectados, t al tiempo, expresado en semanas y c a un coeficiente que mide el numero de contactos por individuo y por unidad de tiempo.
La primera forma de enfocar el problema vendrá dada por la ecuación logística adaptada a este caso:
[math] I'(t)=\frac{c}{N}I(t)(N-I(t)); [/math]
Para seguir el estudio con una aplicación numérica, se fija el número de habitantes N de la ciudad en 500.000. También se sabe que en el instante inicial de la propagación hay 200 personas infectadas y que tras la primera semana se contabilizan 300 nuevos casos.
Estos datos permiten calcular el valor numérico de la constante c mediante el siguiente programa de Octave:
y0=200;
N=500000;
c=linspace(1/100, 99/100, 99);
y=[]; yy=N; yf=N; C=c(1);
for i=1:length(c)
yf=yy;
yy=N*y0/(y0+(N-y0).*exp(-c(i)));
y=[y yy];
if abs(yy-500)<=abs(yf-500)
C=c(i);
end
end
res=sprintf('El valor de c para el que abs(500-y) es minimo (y de valor %d) es %d.', min(abs(y-500)), C);
disp(res);
plot(c,abs(y-500))
xlabel('c');
ylabel('abs(y-500)');
Que se obtiene como el valor de c que minimiza [math] \vert I(primera semana)-500 \vert [/math], que según este programa es 0.92. Resolviendo la ecuación analíticamente el valor de c que se obtiene es 0.917, por lo que se puede considerar que el programa tiene una aproximación aceptable.
Tras este preliminar se dispone de todo lo necesario como para dibujar la curva del número de infectados I(t), para lo que se recurre a aplicar el método de Heun y el método de Runge-Kutta de rango 4.
tf=input('Introduce el numero de semanas a estudiar: ');
h=0.01; c=0.92; N=500000;
t=0:h:tf;
k1=0; k2=0; k3=0; k4=0;
y=[]; yy=200; yyh=200; yh=[];
for i=1:length(t)
k1=c/N*yy*(N-yy);
k2=c/N*(yy+h*k1/2)*(N-yy-h*k1/2);
k2h=c/N*(yy+h*k1)*(N-yy-h*k1);
k3=c/N*(yy+h*k2/2)*(N-yy-h*k2/2);
k4=c/N*(yy+h*k3)*(N-yy-h*k3);
yy=yy+h/6*(k1+2*k2+2*k3+k4);
yyh=yyh+0.5*h*(k1+k2h);
y=[y yy];
yh=[yh yyh];
end
hold on
plot(t,y);
plot(t,yh,'x');
xlabel('Tiempo t (en semanas)');
ylabel('Numero de infectados I(t)');
hold off
res1=sprintf('El numero de infectdos en el tiempo final es %d segun el metodo de Runge-Kutta', yy);
res2=sprintf('El numero de infectdos en el tiempo final es %d segun el metodo de Heun', yyh);
disp(res1); disp(res2);
La epidemia, siguiendo un modelo logístico, crece constantemente y tiende a infectar a la totalidad de los habitantes de la ciudad, que representa una asíntota horizontal para la curva representativa de I(t). Los crecimientos más bruscos del número de infectados van seguidos por un periodo infinito de crecimiento lento (casi nulo) que refleja la dificultad de “cruzarse” con una persona sana a partir de un momento.
Para conocer el número de infectados en un momento concreto, basta con cambiar el programa modificando el tiempo final por el momento a considerar: en la sexta semana hay 45792 infectados.
Con la grafica ya trazada se podrían considerar algunos datos de interés, como por ejemplo a partir de qué momento la epidemia alcanza las 400.000 infecciones (80% de la población), que es 10.05 semanas según:
lim=input('Introduce un numero de infectados: ');
if lim>=500000
while lim>=500000
disp('El numero de infectados debe ser positivo y menor que 500.000.');
lim=input('Introduce un numero de infectados: ');
end
end
h=0.01; c=0.92; N=500000; t0=0; y0=200;
t(1)=0; y(1)=y0;
i=1;
while lim-y(i)>0
y(i+1)=y(i)+h*c/N*y(i)*(N-y(i));
t(i+1)=t(i)+h;
i=i+1;
end
res=sprintf('La epidemia tarda %1.2f semanas en infectar a %d personas', t(end), lim);
disp(res);
plot(t,y)
xlabel('Tiempo (en semanas)');
ylabel('Numero de infectados');
Nota: La condición del bucle if y del bucle while debe ser "lim>=500000|lim<0" pero por alguna razón la wiki no acepta el símbolo y lógico.
Ahora bien, el modelo utilizado implicó muchas simplificaciones iniciales, por lo que no deja de parecer pertinente otro que considere otras situaciones. Se podría pensar en poblaciones tales que si el número de individuos N es elevado entonces la tasa es decreciente, y que si la población en un determinado momento M es demasiado pequeña esta tasa también decrece (por ejemplo, por la dificultad de individuos a los que contagiar los infectados).
De este modo se desarrolla un poco el modelo de la ecuación logística de la forma siguiente: [math] I'(t)=\frac{c}{N}I(t)(N-I(t))(\frac{I(t)}{M}-1); [/math]
y0=[9700 10200 30000];
h=0.01; c=0.92; N=500000; M=10000;
t=0:h:15;
k1=0; k2=0; k3=0; k4=0;
y=[]; Y=zeros(3,length(t));
for k=1:3
yy=y0(k);
y=[y0(k)];
for i=1:length(t)
k1=c/N*yy*(N-yy)*(yy/M-1);
k2=c/N*(yy+h*k1/2)*(N-yy-h*k1/2)*((yy+h*k1/2)/M-1);
k3=c/N*(yy+h*k2/2)*(N-yy-h*k2/2)*((yy+h*k2/2)/M-1);
k4=c/N*(yy+h*k3)*(N-yy-h*k3)*((yy+h*k3)/M-1);
yy=yy+h/6*(k1+2*k2+2*k3+k4);
y=[y yy];
end
y(end)=[];
Y(k,:)=y;
end
hold on
subplot(2,2,1)
plot(t,Y(1,:));
title('Estudio con I(0)=9700');
xlabel('Tiempo (en semanas)');
ylabel('Numero de infectados');
subplot(2,2,2)
plot(t,Y(2,:));
title('Estudio con I(0)=10200');
xlabel('Tiempo (en semanas)');
ylabel('Numero de infectados');
subplot(2,2,3)
plot(t,Y(3,:));
title('Estudio con I(0)=30000');
xlabel('Tiempo (en semanas)');
ylabel('Numero de infectados');
hold off
2 Modelo SIR
Siguiendo la lógica de encontrar un modelo que contemple el mayor número posible de situaciones reales, se considera el modelo SIR.
En esta nueva situación, se divide la población en tres grupos: las personas jamás infectadas (S), las personas infectadas (I), capaces de infectar a otras, y un tercer grupo de personas (R) que engloba a los muertos, las personas que han superado la infección o que están en cuarentena (es decir, las personas infectadas que ya no pueden contagiar a otras personas).
Ahora hacemos suposiciones acerca de los procesos de incubación y de transmisión:
- El número de infectados crece a una tasa proporcional a la cantidad de infectados y a la cantidad de susceptibles: rSI; con r > 0. El número de susceptibles decrece a esta misma tasa, r, llamada tasa de infección.
- La tasa de miembros infectados que pasan a la clase R es proporcional al número de infectados solamente: aI; a > 0. El número de recuperados crece a la misma tasa a, llamada tasa de remoción. Obsérvese que la probabilidad de curarse es una constante, y que podemos asociar 1/a a un tiempo característico de duración de la infección.
- El tiempo de incubación es despreciable, de manera que un susceptible que se infecta se vuelve inmediatamente infeccioso.
Con estas tasas de cambio, la suposición adicional de que la población está bien mezclada, (esto es que todo par de individuos tiene la misma probabilidad de entrar en contacto), y que [math] S(0)=S_0\gt0, I(0)=I_0, R(0)=0, [/math] y por lo tanto [math] N=S_0+I_0 [/math]
podemos escribir el siguiente modelo
[math] \begin{cases} S'=-rSI, \text{ si } t\gt0\\ I'=rSI-aI,\\ R'=I,\\ S(0)=S_0, I(0)=I_0, R(0)=0 \end{cases}[/math]
Observación: [math] S'(t)+I'(t)+R'(t)=0, \forall t\gt0 [/math] Por tanto [math] S(t)+I(t)+R(t) [/math] es constante en todo momento Particularmente en t=0, [math] S(t)+I(t)+R(t)=S(0)+I(0)+R(0)=S_0+I_0=N, \forall t\gt0[/math]
Si nos fijamos en las condiciones iniciales, se aprecia: Si [math] S_0\le \frac{a}{r} [/math], [math] I'_0=rS_0I_0\le r\frac{a}{r}I_0-aI_0=0 [/math] La derivada en el instante inicial es negativa, por lo que la epidemia tiende a desaparecer en vez de a expandirse.
Vamos a proceder a representar nuestro modelo generalizado, aplicando los metodos de Euler y Runge-Kutta, para distintos valores de [math] S_0 ,I_0 [/math] y h como queda reflejado en las graficas siguientes:
S0=[159985 140000 159999]; I0=[15 20000 1]; %vectores de valores iniciales
hold on
for i=1:3
w0=[S0(i); I0(i); 0]; %vector inicial del sistema
h=0.01; %paso
t0=0; tN=20; %intervalo de tiempo a estudiar
N=(tN-t0)/h;
t=t0:h:tN;
r=0.0000218; a=0.341; %parametros del modelo
w=zeros(3,N+1); y=zeros(3,N+1);
w(:,1)=w0; y(:,1)=w0;
ww=w0; yy=w0; A=[0 0 0; 0 -a 0; 0 a 0]; %matriz del sistema
for n=1:N
F=[-r*ww(1)*ww(2); r*ww(1)*ww(2); 0];
k1=A*yy+F;
k2=A*(yy+h*k1/2)+F;
k3=A*(yy+h*k2/2)+F;
k4=A*(yy+h*k3)+F;
yy=yy+h/6*(k1+2*k2+2*k3+k4);
ww=ww+h*(A*ww+F);
w(:,n+1)=ww; %vector de calculo con el método de Euler
y(:,n+1)=yy; %vector de calculo con el método de R-K
end
subplot(2,3,i);
plot(t,w,'Linewidth',4); %graficas calculadas por el método de Euler
titulo=sprintf('Euler S0=%d, I0=%d', S0(i), I0(i));
legend('S(t)','I(t)', 'R(t)');
title(titulo);
xlabel('Tiempo t (en semanas)');
subplot(2,3,i+3);
plot(t,w,'Linewidth',4); %graficas calculadas por el método de R-K
titulo=sprintf('R-K S0=%d, I0=%d', S0(i), I0(i));
legend('S(t)','I(t)', 'R(t)');
title(titulo);
xlabel('Tiempo t (en semanas)');
end
Podemos ver como al reducir el valor de nuestra h, vamos reduciendo el error cometido y nos acercamos más a la solución real.
Podemos observar en los casos estudiados con anterioridad, en los cuales varía el valor inicial, como aquellas personas susceptibles de ser infectadas, representadas por la letra S, y las infectadas, por la letra I, evolucionan de modo inverso, es decir, mientras los susceptibles van decreciendo en función del tiempo, los infectados y removidos van creciendo, quedando demostrado que las variables utilizadas, se ajustan a los valores numéricos. Así, si ahora tomamos una porción de tiempo determinado, se observan en las gráficas obtenidas, como ocurren variaciones, siendo a causa de que el número de personas susceptibles a ser infectadas empieza a decrecer o debido a que se llega a un momento en el que las variables tienden a igualarse.
Empezaremos por representar gráficamente la ecuación [math] I+S-\rho log(S)=k [/math] que viene de la ecuación diferencial dada en el enunciado.
S0=[159985 140000 159999]; I0=[15 20000 1];
colores=['r' 'b' 'k']; Leyenda=[];
hold on
for i=1:3
y0=[S0(i); I0(i); 0];
h=0.01;
t0=0; tN=20;
N=(tN-t0)/h;
t=t0:h:tN;
r=0.0000218; a=0.341;
y=zeros(3,N+1);
y(:,1)=y0;
yy=y0; A=[0 0 0; 0 -a 0; 0 a 0];
for n=1:N
F=[-r*yy(1)*yy(2); r*yy(1)*yy(2); 0];
k1=A*yy+F;
k2=A*(yy+h*k1/2)+F;
k3=A*(yy+h*k2/2)+F;
k4=A*(yy+h*k3)+F;
yy=yy+h/6*(k1+2*k2+2*k3+k4);
y(:,n+1)=yy;
end
plot(y(1,:),y(2,:),colores(i),'Linewidth',6);
k=y(2,:)+y(1,:)-(0.341/0.0000218)*log(y(1,:));
leyenda=sprintf('k=%1.3f', k(1));
Leyenda=[Leyenda leyenda];
title('Trayectorias I+S-rho*log(S)=k');
xlabel('Numero de personas sanas S(t)');
ylabel('Numero de personas infectadas I(t)');
end
legend(Leyenda(1:12),Leyenda(13:24),Leyenda(25:36));
Para empezar el intervalo de valores en el que se encuentra la gráfica, es aceptable. Ya que los vértices del triángulo serían los casos ideales en los que; para el total de la población N, que estaría toda sana, el número de infectados sería 0; de ahí que sea (N,0). Para el caso contrario se deduce lo mismo para el punto (0,N).
Evidentemente, el resultado que obtenemos en la gráfica se encuentra dentro del intervalo indicado. Otro aspecto que cabe destacar de la gráfica es que las soluciones gráficas son "cuasiparalelas" ya que la diferencia entre ellas se ve reflejada en el valor inicial que desplaza las gráficas siguientes respecto de la inicial. Como la diferencia entre los primeros valores iniciales hay una diferencia menor que la de ambas con respecto a la tercera, en la gráfica tendría que aparecer representado. Se puede ver a simple vista que esto ocurre ya que las dos primeras se ven casi iguales y la tercera está más diferenciada.
Ahora procederemos a calcular en qué semana habrá un número máximo de infectados.
S=y(2,:); maximo=0;
tiempo=0; m=0;
for l=2:length(S)
m=m+h;
if S(l)>S(l-1)
tiempo=m;
maximo=S(l);
end
end
tiempo
maximo
Para el cálculo del máximo, hemos decidido que la mejor manera de llevarlo a cabo, era representar gráficamente la función que nos indica el enunciado [math] I(maximo)=N-\rho+\rho log \frac{\rho}{S_0} [/math]
Una vez obtenemos la gráfica, con el comando max obtenemos el valor numérico de dicho máximo y como se puede observar coinciden.
Max=10817 en 4.62 semanas.
