Revisión del 16:08 6 mar 2015

1 Introducción y consideraciones iniciales

Se va a proceder al estudio del comportamiento temporal de una enfermedad infecciosa sin extensión espacial.
El crecimiento de una epidemia podrá asimilarse al modelo logístico.

1.1 Hipótesis iniciales

1. La población es un número fijo y cada miembro de la población es susceptible a la enfermedad.
2. La duración de la enfermedad es larga, de manera que no se cura durante el periodo de estudio.
3. Todos los individuos infectados son contagiosos y circulan libremente entre la población.
4. Cada contacto de una persona infectada con una persona no infectada redunda en la transmisión de la enfermedad.

1.2 Valores iniciales

• La ciudad de estudio posee 500.000 habitantes.
• Al inicio de la primera semana de registro se contabilizaron 200 casos.
• Durante la primera semana aparecieron 300 nuevos casos.

1.3 Variables

• t: Periodo de tiempo de análisis (en semanas).
• N: Población.
• c: Número de contactos con otros individuos de cada persona infectada por unidad de tiempo.

2 Modelo EDO de primer orden

Se llamará I(t) al número de infectados al tiempo t.
El problema vendrá dada por la ecuación logística: [math] I'(t)=\frac{c}{N} \cdot I(t) \cdot (N-I(t)) [/math]

2.1 Número de contactos de una persona infectada

Utilizaremos un método numérico para calcular mediante Euler el número de contactos de una persona infectada con otros individuos por semanas. Elegiremos el valor c aquel que minimice los valores calculados, mediante la siguiente fórmula:

[math] \vert I (primera semana) − 500 \vert [/math]

Fig 1: Relación entre el valorc y el error cometido.

h=0.01; 
c0=0.01; 
cN=0.99;
c=c0:h:cN;
k=length(c);
y0=200;
y(1)=y0;
N=500000;
for i=1:99;
    c(i)=i/100;
    m(i)=N*y0/(y0+(N-y0).*exp(-c(i)));
    z(i)=abs(500-m(i));
end
plot(c,z)

2.2 Evolución de la epidemia

Deseamos conocer la evolución del número de infectados a lo largo de las semanas, para ello emplearemos los métodos de Heun y Runge-Kutta de orden 4.

Fig 2: Gráfica de infectados según el método de Heun.

N=500000;
h=0.01;
I=[];
I0=200;
I(1)=I0;
c=0.92;
n=1;
t0=0;
while I<400000
    % Método de Heun
K1=I(n)*(c/N)*(N-I(n));
K2=(I(n)+K1*h)*(c/N)*(N-I(n)+K1*h);
I(n+1)=I(n)+(h/2)*(K1+K2);
n = n+1;
end
tNI=h*(length(I)-1);
tI=t0:h:tNI; 
plot(tI,I)
valor6semanasheun= I(6/h)
tiempo400000infectados=round(n/100)

Fig 3: Gráfica de infectados según el método de Runge-Kutta.

N=500000;
h=0.01;
s0=200;
s(1)=s0;
c=0.92;
k=1;
t0=0;
% Método Runge-Kutta
while s<400000
C1=s(k)*(c/N)*(N-s(k));
C2=(s(k)+1/2*C1*h)*(c/N)*(N-(s(k)+1/2*C1*h));
C3=(s(k)+1/2*C2*h)*(c/N)*(N-(s(k)+1/2*C2*h));
C4=(s(k)+C3*h)*(c/N)*(N-(s(k)+h*C3));
s(k+1)=s(k)+(h/6)*(C1+2*C2+2*C3+C4);
k=k+1;
end
tNS=h*(length(s)-1);
tS=t0:h:tNS;
plot(tS,s,'g')
Valor6semanasrungekutta= s(6/h)
tiempo400000infectados=round(k/100)

Comprobamos que en las gráficas el número de infectados es siempre creciente.

No podemos definir un intervalo de máxima infección ya que las gráficas crecen de forma exponencial, es decir el número de infectados seguirá incrementándose hasta que el total de la población haya sido infectada.

Además de las gráficas, el código nos devuelve el número de infectados para 6 semanas y el tiempo que se tarda en llegar a 400.000 infectados.

• Método de Heun: A las 6 semanas hay 45.065 infectados. Se llega a los 400.000 infectados a las 10 semanas.
  

• Método de Runge-Kutta: A las 6 semanas hay 45.032 infectados. Se llega a los 400.000 infectados a las 10 semanas.

Fig 4: Evolución de la infección para los 3 valores iniciales.

N= 500000; % población total
M= 10000; % umbral de población que asigna tasas de infección decrecientes ara valores inferiores
c=0.92; % contagio que optimiza el error(apartado 1)
ti=0; tf=4*3; % calcularemos la gráfica para un período de 3 meses de propagación (4 semanas por 12 meses)
h=0.01;
t= ti:h:tf;
k=(tf-ti)/h;
I1=zeros(1,k+1);
I10=9700; 
I1(1)=I10;

for i=1:k
    K1 = (c*I1(i))/N*(N-I1(i))*((I1(i)/M)-1);
    K2 = (c*(I1(i)+K1*h))/N*(N-(I1(i)+K1*h))*(((I1(i)+K1*h)/M)-1);
    I1(i+1) = I1(i) + (h/2)*(K1+K2);
end
I2=zeros(1,k+1);
I20=10200; 
I2(1)=I20;

for j=1:k
    L1 = (c*I2(j))/N*(N-I2(j))*((I2(j)/M)-1);
    L2 = (c*(I2(j)+L1*h))/N*(N-(I2(j)+L1*h))*(((I2(j)+L1*h)/M)-1);
    I2(j+1) = I2(j) + (h/2)*(L1+L2);
end
I3=zeros(1,k+1);
I30=30000;
I3(1)=I30;

for l=1:k
    P1 = (c*I3(l))/N*(N-I3(l))*((I3(l)/M)-1);
    P2 = (c*(I3(l)+K1*h))/N*(N-(I3(l)+P1*h))*(((I3(l)+P1*h)/M)-1);
    I3(l+1) = I3(l) + (h/2)*(P1+P2);
end
hold on
plot(t,I1)
plot(t,I2,'m')
plot(t,I3,'g')
legend('Io=9700','Io=10200','Io=30000','location','best')
hold off

Fig 5: Infectados iniciales Io=9.700.

Los infectados solo tienden a desaparecer en el caso de [math] I_0=9.700 [/math], pues [math] I_0\ltM [/math]. En los otros dos casos la infección crece por ser [math] I_0\gtM [/math].

Fig 6: Infectados iniciales Io=10.200.

La evolución de los infectados en el caso en los que Io sea igual a 10.200, alcanza su máximo aproximadamente en la semana 4. y en el caso en que [math] I_0=30.000 [/math], alcanza su máximo aproximadamente en la
semana 1.

Fig 7: Infectados iniciales Io=30.000.

En el caso de [math] I_0=9.700 [/math], al ser una evolución decreciente el máximo se sitúa al inicio del estudio.

3 Modelo SIR

Vamos a considerar la utilización del modelo SIR para poder abarcar un estudio más general.
Para ello supondremos una población total constante de N individuos, así como la introducción de un grupo de infectados que es introducido en la población.
Plantearemos así la descripción temporal del grupo de individuos infectados.

3.1 Clases de población

Impondremos en el estudio que los individuos que se recuperen de la enfermedad, se vuelven inmunes a ella, resultando tres clases de población:

• Individuos susceptibles de infectarse: S
  
• Individuos ya infectados con capacidad de infectar a individuos sanos: I
  
• Grupo de individuos recuperados, inmunes, aislados y muertos: R
  

[math] S \rightarrow I \rightarrow R [/math]

Llamaremos S(t), I(t) y R(t) a la cantidad de individuos susceptibles de infectarse, ya infectados y del grupo de recuperados, inmunes, aislados y muertos respectivamente; al cabo de t unidades de tiempo.

3.2 Suposiciones sobre el desarrollo de la enfermedad

• El número de infectados crece proporcionalmente al número de infectados y susceptibles, denominaremos a esta tasa de infección r. El número de individuos susceptibles de infectarse decrece con esta misma tasa r.
• La tasa de miembros infectados que pasan a la clase de removidos es proporcional al número de infectados solamente. El número de individuos recuperados crece con una tasa a, que denominaremos tasa de remoción.
• El tiempo de incubación es despreciable, es decir un individuo que siendo susceptible de padecer la infección, se infecta, automáticamente se vuelve individuo infeccioso.

3.3 Demostraciones realizadas

El sistema que gobierna el modelo SIR planteado será:

[math] \begin{cases} S'=-rSI, \text{ para } t\gt0\\ I'=rSI-aI,\\ R'=I,\\ S(0)=S_0, I(0)=I_0, R(0)=0 \end{cases} [/math]

Vamos a demostrar que [math] S(t)+I(t)+R(t)=S_0 +I_0 =N [/math], para un tiempo mayor o igual a 0.

A la izquierda de la igualdad tenemos la suma de tres funciones que dependen del tiempo, mientras que a la derecha tenemos la suma de dos valores constantes.

Como consecuencia de esto, si derivamos en ambos lados respecto del tiempo, tendremos que la suma de las derivadas de las funciones tendrá que ser nula:

[math] S’(t)+I’(t)+R’(t)=0 : (Ec. 1) [/math]

Según las ecuaciones planteadas en el sistema indicado:

[math] S’(t)= -r \cdot S(t) \cdot I(t) : (Ec.2) [/math]

[math] I’(t)=r \cdot S(t) \cdot I(t)-a \cdot I(t) : (Ec.3) [/math]

[math] R’(t)=a \cdot I(t) : (Ec.4) [/math]

Sustituyendo en la Ec. 1:

[math] S’(t)+I’(t)+R’(t)= -r \cdot S(t) \cdot I(t)+r \cdot S(t) \cdot I(t)-a \cdot I(t)+a \cdot I(t)=0 [/math]

Quedando demostrada la igualdad.

A continuación vamos a demostrar como cuando [math] S_0 \le \displaystyle\frac{a}{r} [/math], no hay epidemia ya que [math] I(t) \le I_0 [/math]. Es decir, el número de infectados no aumenta y tienden a desaparecer de forma breve.

Empezando por el lado derecho de la inecuación:

Si [math] I(t) \le I_0 [/math], derivando respecto del tiempo y sabiendo que [math] I(0)=I_0 [/math] constante, vemos que [math] I'(t) \le 0 [/math]

De la ec. 3 tenemos:

[math] r \cdot S(t) \cdot I(t)-a \cdot I(t) \le 0 \longrightarrow r \cdot S(t) \cdot I(t) \le a \cdot I(t) \longrightarrow S(t) \le \displaystyle\frac{a}{r} [/math]

Al cumplirse esta inecuación para todo t, queda demostrada la desigualdad.

Por último demostraremos como si [math] S_0 \gt \displaystyle\frac{a}{r} [/math] hay epidemia ya que el número de infectados aumentan, es decir [math] I(t) \gt I_0 [/math].

Operando de forma similar al caso anterior tenemos:

[math] r \cdot S(t) \cdot I(t)-a \cdot I(t) \gt 0 \longrightarrow r \cdot S(t) \cdot I(t) \gt a \cdot I(t) \longrightarrow S(t) \gt \displaystyle\frac{a}{r} [/math]

3.4 Cálculo numérico de la epidemia

Plantearemos un cálculo numérico mediante los métodos de Euler, Euler implícito y Runge-Kutta de orden 4 para tres casos:

• Caso 1: [math] S_0 = 159985, I_0 = 15, r = 0.0000218, a = 0.341, I = [0,4] y [0,20] [/math]
• Caso 2: [math] S_0 = 140000, I_0 = 20000, r = 0.0000218, a = 0.341, I = [0,4] y [0,20] [/math]
• Caso 3: [math] S_0 =159999, I_0 =1, r=0.0000218, a=0.341, I=[0,4] y [0,20] [/math]

t0=0;h= 0.1;tf=4; % el tiempo total es conocido, calculamos gráficas para los intervalos [0,4] y [0,20] y el paso h también cambia
S0=159985; I0=15; % datos de la población (los vamos cambiando para cada gráfico)
r=0.0000218; a=0.341; R0=0;
t=t0:h:tf;
k=((tf-t0)/h);
Sa=[]; %el subíndice a se relaciona con el método de euler
Ia=[];
Ra=[];
Sb=[]; % el subíndice b con euler modificado
Ib=[];
Rb=[];
Sc=[]; % el c con runge-kutta
Ic=[];
Rc=[];
Sa(1)=S0;
Ia(1)=I0;
Ra(1)=R0;
Sb(1)=S0;
Ib(1)=I0;
Rb(1)=R0;
Sc(1)=S0;
Ic(1)=I0;
Rc(1)=R0;
% euler
for i=1:k
    Ia(i+1)=Ia(i)+h*(r*Sa(i)*Ia(i)-a*Ia(i));
    Sa(i+1)=Sa(i)+h*(-r*Sa(i)*Ia(i));
    Ra(i+1)=Ra(i)+h*(a*Ia(i));
end
% euler modificado
for j=1:k
    KI=r*Sb(j)*Ib(j)-a*Ib(j);
    Ib(j+1)=Ib(j)+h*(r*Sb(j)*(Ib(j)+h*0.5*KI)-a*(Ib(j)+h*0.5*KI));
    KS=-r*Sb(j)*Ib(j);
    Sb(j+1)=Sb(j)+h*(-r*(Sb(j)+h*0.5*KS)*Ib(j));
    KR=a*Rb(j);
    Rb(j+1)=Rb(j)+h*(a*(Ib(j)));
end
% runge-kutta
for n=1:k
 k1S = -r*Sc(n)*Ic(n); 
 k1I = r*Sc(n)*Ic(n)-a*Ic(n);
 k1R = a*Ic(n);
 k2S = -r*(Sc(n) + 1/2*k1S*h)*(Ic(n) + 1/2*k1I*h);
 k2I = r*(Sc(n) + 1/2*k1S*h)*(Ic(n) + 1/2*k1I*h)-a*(Ic(n) + 1/2*k1I*h);
 k2R = a*(Ic(n) + 1/2*k1I*h);
 k3S = -r*(Sc(n) + 1/2*k2S*h)*(Ic(n) + 1/2*k2I*h);
 k3I = r*(Sc(n) + 1/2*k2S*h)*(Ic(n) + 1/2*k2I*h)-a*(Ic(n) + 1/2*k2I*h);
 k3R = a*(Ic(n) + 1/2*k2I*h);
 k4S = r*(Sc(n) + k3S*h)*(Ic(n) + k3I*h);
 k4I = r*(Sc(n) + k3S*h)*(Ic(n) + k2I*h)-a*(Ic(n) + 1/2*k2I*h);
 k4R = a*(Ic(n) + k3I*h);
Sc(n+1) = Sc(n) + (h/6)*(k1S+2*k2S+2*k3S+k4S); 
Ic(n+1) = Ic(n) + (h/6)*(k1I+2*k2I+2*k3I+k4I); 
Rc(n+1) = Rc(n) + (h/6)*(k1R+2*k2R+2*k3R+k4R); 
end
hold on
plot(t,Ia,'--g',t,Sa,'--m',t,Ra,'--b')
plot(t,Ib,':g',t,Sb,':m',t,Rb,':b')
plot(t,Ic,'-g',t,Sc,'-m',t,Rc,'-b')
legend('I Euler','S Euler','R Euler','I Euler modificado','S Euler modificado','R Eler modificado','I Runge-Kutta','S Runge-Kutta','R Runge-Kutta','Location','best')
hold off

EDG17_2015_apartado31h1

EDG17_2015_apartado31h1b

EDG17_2015_apartado31h2

EDG17_2015_apartado31h2b

EDG17_2015_apartado32h1

EDG17_2015_apartado32h1b

EDG17_2015_apartado32h2

EDG17_2015_apartado32h2b

EDG17_2015_apartado33h1

EDG17_2015_apartado33h1b

EDG17_2015_apartado33h2

EDG17_2015_apartado33h2b

3.5 Interpretación de la epidemia

Consideraremos que la epidemia ha finalizado cuando el número de infectados deja de crecer, es decir, cuando la función I(t) alcanza su máximo absoluto.

Podemos ver claramente que todos los casos corresponden a modelos de epidemias, ya que siempre se cumple que los distintos valores de [math] S_0 [/math] son superiores al coeficiente de la tasa de remoción relativa [math](\rho = \frac{a}{r} ) [/math]

NOTA: En este apartado hemos representado cada gráfica solamente con dos pasos en vez de tres, ya que como se ha podido ver la diferencia de precisión es imperceptible, de manera que para no saturar con demasiadas gráficas hemos optado por obviar los casos en los que h=0.01.

En caso de que quisiésemos conocer estas gráficas bastaría con copiar el código de las anteriores sustituyendo el valor de h.

En todos los casos la diferencia de tamaños de los intervalos escogidos nos permiten comprobar mejor cómo es el comportamiento de las diferentes variables ya que en algunos casos, hay intervalos que no tienen la longitud suficiente para recoger toda la información. Especialmente un intervalo de mayor tamaño nos permite saber a qué valores tienden asintóticamente el número de ciudadanos infectados, susceptibles y removidos al cabo de un largo periodo de tiempo.

Interpretaremos por tanto los 3 casos propuestos en el apartado anterior:

Caso 1:

• De acuerdo con los datos analíticos, podríamos estimar que la epidemia finaliza cuando el número de infectados comienza a decrecer, ya que eso nos indica que cada vez menos personas contraen la enfermedad, y como podemos ver en las gráficas, el número de personas susceptibles de contagiarse ha empezado a decrecer desde que se inició la epidemia.
• Decimos que se cumple que la duración media de la epidemia sea de [math] \frac{1}{a} [/math] (tres semanas aproximadamente) pues vemos claramente que las curvas de infectados alcanzan su máximo en torno a la cuarta semana desde que comienza la epidemia, aunque podríamos considerarlo como una estimación.
• El número de infectados no es inferior al inicial hasta la semana 20.
• De acuerdo con la información del gráfico, la tasa de reproducción de la infección debería ser inferior.

Caso 2:

• A pesar de que el número de infectados iniciales es considerablemente mayor que en el caso anterior, vemos cómo la gráfica sigue un comportamiento epidémico análogo al mostrado en el caso anterior, aunque con ligeras diferencias.
• En este caso no podemos considerar la duración media de la epidemia como [math] \frac{1}{a} [/math], ya que como podemos ver en las gráficas el número de infectados empieza a decrecer aproximadamente al cabo de 9 días.
• También aquí vemos que la tasa de reproducción de la infección [math] (R_0=\frac{r}{a} \cdot S_0)[/math] debería ser menor que el valor obtenido analíticamente, aunque sin duda el valor obtenido previo a la gráfica nos da una idea previa de cómo va a ser la diferencia de tamaño entre los habitantes susceptibles iniciales y los removidos.

Caso 3:

• El comportamiento es muy similar al caso 1, solo que con un ligero desplazamiento a la derecha. Podríamos considerar también que la duración media de la epidemia es de 3 semanas, aunque se ve claramente que no es hasta la quinta semana cuando empieza a disminuir el número de infectados.

Como conclusión final vamos a destacar la relación entre el número inicial de habitantes y el valor al que tiende la función R(t) al cabo del tiempo. Con el transcurso de las semanas vemos en todas las gráficas que el número de personas removidas tiende a 250000, es decir, a N/2. Por tanto la mitad de la población ha sido afectada por la epidemia.

3.6 Trayectorias