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.

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
min(abs(y-500))
C
plot(c,abs(y-500))

Que se obtiene como el valor de c que minimiza [math] abs{I(primera semana)-500} [/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.

h=0.001; c=0.92; N=500000;
t=0:h:20;
k1=0; k2=0; k3=0; k4=0;
y=[]; yy=200;
for i=1:length(t)
  k1=c/N*yy*(N-yy);
  k2=c/N*(yy+h*k1/2)*(N-yy-h*k1/2);
  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);
  y=[y yy];
end
plot(t,y);
yy

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.

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.4 semanas según:

h=0.1; c=0.92; N=500000; t0=0; y0=200;
t(1)=0; y(1)=y0;
i=1;
lim=400000;
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
t(end)
plot(t,y)

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))(I/M-1); [/math]

h=0.001; c=0.92; N=500000; M=10000;
t=0:h:15;
k1=0; k2=0; k3=0; k4=0;
y=[]; yy=9700;
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
plot(t,y);