Diferencia entre revisiones de «Modelos epidemiológicos (Grupo 3A)»

Revisión actual del 08:48 12 may 2016

En el desarrollo de una epidemia, se distinguen dos tipos de individuos. Los que ya han contraído la enfermedad (infectados), que llamaremos [math]I(t)[/math] ; y los que son susceptibles de contraerla (sanos), a los que llamaremos [math]S(t)[/math]. Donde [math]t[/math] es la variable del tiempo. Se presentan dos hipótesis para realizar este estudio:

1. La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del número de personas infectadas.

2. La tasa de individuos que pasan de ser susceptibles a contraer la enfermedad a estar infectados es proporcional a la interacción entre el número de individuos en ambas clases.

A través del siguiente sistema de ecuaciones diferenciales en función del tiempo observamos las variaciones de ambas poblaciones:

[math]\begin{cases}\frac{dS}{dt}=-aSI\\\frac{dI}{dt}=aSI-bI-cI\end{cases}[/math]

Donde [math]a[/math], [math]b[/math] y [math]c[/math] son parámetros.

• [math]\frac{dS}{dt}[/math] : Es la variación de la población susceptible de contraer la enfermedad.

• [math]\frac{dI}{dt}[/math] : Es la variación de la población infectada.

• [math]aSI[/math] : Es la interacción entre ambas poblaciones que dependen de [math]a[/math].

• [math]bI[/math] : Indica los individuos que fallecen. Por lo tanto [math]b[/math] podría considerarse como la tasa de curación.

• [math]cI[/math] : Indica los individuos que se han curado. Por lo tanto [math]c[/math] podría considerarse como la tasa de mortalidad.

Para interpretar [math]a[/math] observamos el sistema y la segunda hipótesis. Como los individuos que pasan de ser susceptibles a estar infectados son proporcionales a su interacción [math](SI)[/math], [math]a[/math] podría considerarse la tasa de contagio en la población susceptible.

1 Resolución del problema con una sola ecuación diferencial

Primero vamamos a estudiar el modelo epidemiológico de una población en la que incialmente las 2000 personas que conforman la población se encuentran infectadas. Para hacer el seguimiento de la enfermedad contamos tanto con la constante de proporcionalidad de personas que se curan (b=0.3) y la constante de proporcionalidad del número de personas que fallecen (c=0.01). De esta forma nos quedará el siguiente PVI:

[math] \begin{cases}   \frac{dI}{dt} = -0.31  I \\I_{0}  = 2000\end{cases}  [/math]

A continuación vamos a proceder a resolverlo por el método de Euler y trapecio asignandole un tamaño de paso de h=0.1 y veremos cuanto tiempo tardará en reducirse el número de infectados a la cuarta parte.

[math]\Rightarrow[/math] Método de Euler

Aproximación mediante Euler

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
b=0.3;
c=0.01;

t(1)=t0;
ye(1)=y0; %Vector solución Euler

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500

%EULER

while ye(i)>500
    
    ye(i+1)=ye(i)+h*(-b*ye(i)-c*ye(i)); 
    t(i+1)=t(i)+h;
    i=i+1;
end

[t',ye'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,ye)
legend('Población infectada(Método de Euler)','Location','best');

[math]\Rightarrow[/math] Método del trapecio

Aproximación mediante el método del trapecio

clear all
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0; %Vector solución trapecio

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
while z(i)>500
     z(i+1)=(z(i)*(1-(h/2)*(b+c)))/(1+h/2*(b+c));%trapecio
     t(i+1)=t(i)+h;
     i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada(Método del trapecio)','Location','best');

Con todo esto podemos determinar que el número de infectados se reducirá a la cuarta parte en 4.5 días, esto ocurrirá o bien por el fallecimiento o por la curación de dicha enfermedad. Cabe señalar que el número de curados será mayor que el de defunciones, esto se puede ver a simple vista viendo que b>c. De cara a nuestro estudio de esta enfermedad vamos a ver como afectaría al desarrollo de la misma si introducimos 100 sujetos sanos en la población, tomando una constante de interacción entre personas infectadas y sanas de a=0.003. Procederemos igual que antes a resolverlo por el método de Euler y trapecio con un tamaño de paso de 0.1. Primero vamos a definir nuestro nuevo PVI:

[math] \begin{cases}   \frac{dI}{dt} = -0.01  I \\I_{0}  = 2000\end{cases}  [/math]

[math]\Rightarrow[/math] Método de Euler

Aproximación mediante Euler para S=100

%Euler
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
y0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
y(1)=y0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de y(i) hasta que y(i)>500
while y(i)>500
    
    y(i+1)=y(i)+h*(a*S*y(i)-b*y(i)-c*y(i));
    t(i+1)=t(i)+h;
    i=i+1;
end

[t',y'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,y)
legend('Población infectada para S=100(Método de Euler)','Location','best');

[math]\Rightarrow[/math] Método del trapecio

Aproximación mediante el método del trapecio para S=100

%Trapecio
clear all
%Datos
%S: Es la población susceptible de ser infectados
%I: Es la población de individuos infectada
%a: Es la constante de proporcionalidad entre los infectados por cada
%interacción
%b: Es la tasa de fallecimientos de las personas infectadas
%c: Es la tasa de las personas que se curan de las infectadas

%Damos valor al tamaño de paso
h=0.1;

%Definimos las constantes y las condiciones iniciales
t0=0;
z0=2000;
S=100;
a=0.003;
b=0.3;
c=0.01;

t(1)=t0;
z(1)=z0;

i=1; % Lo usaremos a continuación en el bucle

%Con el bucle while para calcularemos valores de z(i) hasta que z(i)>500

%TRAPECIO
  while z(i)>500
      z(i+1)=(z(i)*(1-(h/2)*(b+c-S*a)))/(1+h/2*(b+c-S*a));
      t(i+1)=t(i)+h;
      i=i+1;
end

[t',z'];

disp('Tiempo en el que el número de infectados se reduce a la cuarta parte')
disp(t(end))
plot(t,z)
legend('Población infectada para S=100(Método del trapecio)','Location','best');

Ahora con estos nuevos datos vemos que para que se reduzca a la cuarta parte el número de infectados tendrán que pasar 138.7 días, algo lógico al fijarnos en la ED, ambas son ecuaciones decrecientes pero para el caso de S=100 la pendiente será mucho menor. Si continuamos haciendo pruebas vemos para valores superiores de s, por ejemplo, S=200 el matlab no nos da ninguna respuesta. Esto ocurre porque al aumentar S la ecuación pasa a ser creciente, por lo que lógicamente si intentas buscar un valor del número de infectados menor al inicial no lo podrás encontrar.

2 Resolución del problema completo

Vamos a resolver el problema completo para ver como evolucionarían las poblaciones en un periodo de 40 días. Primero supongamos una población infectada inicial de 20 individuos y una población en riesgo de contagio de 800. Después supondremos el caso de 40 individuos enfermos y 10000 en riesgo de contagio. Para resolverlo usaremos el método de Euler y más tarde el método de Runge-Kutta de orden cuatro con distintas discretizaciones para ver la influencia de estas.

[math]\Rightarrow[/math] Método de Euler

%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:'); 
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;          
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
    y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h
legend('Población sana','Location','best','Población enferma','Location','best');

[math]\Rightarrow[/math] Método de Runge-Kutta

%Runge-Kutta
clear all

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
a=0.003;
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
     K1S=-a*y(1,i)*y(2,i);
     K2S=-a*(y(1,i)+(1/2)*K1S*h)*(y(2,i)+(1/2)*K1S*h);
     K3S=-a*(y(1,i)+(1/2)*K2S*h)*(y(2,i)+(1/2)*K2S*h);
     K4S=-a*(y(1,i)+K3S*h)*(y(2,i)+K3S*h);
     y(1,i+1)=y(1,i)+(h/6)*(K1S+2*K2S+2*K3S+K4S);
     K1I=a*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
     K2I=a*(y(1,i)+(1/2)*K1I*h)*(y(2,i)+(1/2)*K1I*h)-b*(y(2,i)+(1/2)*K1I*h)-c*(y(2,i)+(1/2)*K1I*h);
     K3I=a*(y(1,i)+(1/2)*K2I*h)*(y(2,i)+(1/2)*K2I*h)-b*(y(2,i)+(1/2)*K2I*h)-c*(y(2,i)+(1/2)*K2I*h);
     K4I=a*(y(1,i)+K3I*h)*(y(2,i)+K3I*h)-b*(y(2,i)+K3I*h)-c*(y(2,i)+K3I*h);
     y(2,i+1)=y(2,i)+(h/6)*(K1I+2*K2I+2*K3I+K4I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Gráfico
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
legend('Poblacion sana','Poblacion enferma','Location','best'); 
hold off

% vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

La dificultad que hay en el uso del método del trapecio o cualquier otro método implícito está en que las incógnitas de nuestro sistema [math](S,I)[/math] quedan implícitas en ambas ecuaciones y necesitan que se despejen. Al comenzar a realizar cálculos nos damos cuenta que el despeje es imposible, ya que las dos variables incógnita aparecen multiplicándose entre sí, quedando una variable dependiendo de la otra siempre.

2.1 Situación inicial de [math](S _{0} , I_{0} )=(800,20)[/math]

En estos programas, sustituyendo las [math]h[/math] para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para [math]h=0.1[/math]

Método de Euler

Método de Runge-Kutta

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 prácticamente toda la población ha fallecido.

• Método de Euler: nuestro programa nos indica que en este caso el momento en el que los enfermos son máximos se produce al tercer día (2.8) y el número de enfermos en ese momento es 517.
• Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce al tercer día (2.6) y el número de enfermos en ese momento es 650.

-Para [math]h=0.01[/math]

Método de Euler

Método de Runge-Kutta

En el gráfico apenas se encuentra diferencia, lo cual es lógico porque es el mismo problema, pero si sobrepusiéramos las gráficas, lograríamos ver como se va ajustando a la realidad en función del aumento de h.

• Método de Euler: en este caso el momento en el que los enfermos son máximos es también en el tercer día (2.7) y el número de enfermos es 506.
• Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce al tercer día (2.69) y el número de enfermos en ese momento es 517.

-Para [math]h=0.001[/math]

Método de Euler

Método de Runge-Kutta

• Método de Euler: en este caso el momento en el que los enfermos son máximos se produce el tercer día (2.696) con un valor de 505 enfermos.
• Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce el tercer día (2.694) con un valor de 506 enfermos.

-Para [math]h=0.0001[/math]

Método de Euler

Método de Runge-Kutta

• Método de Euler: al igual que antes, se repite el momento en el que los enfermos son máximos en el tercer día (2.695) con un valor de 505 enfermos.
• Método de Runge-Kutta: se repite el momento en el que los enfermos son máximos en el tercer día (2.6948) con un valor de 505 enfermos.

Observamos en las distintas gráficas obtenidas, que según vamos empleando un tamaño de paso [math]h[/math] más pequeño, nos vamos acercando cada vez más a la población real que pretendemos reflejar. Nos damos cuenta de que las gráficas tienen bastante lógica, pues a la vez que desciende el número de población sana va aumentando la población infectada. También observamos que para el método de Runge-Kutta cuando utilizamos tamaño de paso [math]h=0.1[/math] pierde precisión en relación a cuando utilizamos tamaños de paso más pequeños.

2.2 Situación inicial de [math](S _{0} , I_{0} )=(10000,40)[/math]

En estos programas, sustituyendo las [math]h[/math] para 0.1, 0.01, 0.001 y 0.0001 obtendremos los siguientes resultados:

-Para [math]h=0.1[/math]

Método de Euler

Método de Runge-Kutta

• Método de Euler: el momento en el que los enfermos son máximos se produce el primer día (0.5) con un valor de 12600 enfermos.
• Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce el primer día (0.4) con un valor de 2.6e74 enfermos.

-Para [math]h=0.01[/math]

Método de Euler

Método de Runge-Kutta

• Método de Euler: el momento en el que los enfermos son máximos se produce el primer día (0.35) con un valor de 9521 enfermos.
• Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce el primer día (0.31) con un valor de 13550 enfermos.

-Para [math]h=0.001[/math]

Método de Euler

Método de Runge-Kutta

• Método de Euler: el momento en el que los enfermos son máximos se produce el primer día (0.342) con un valor de 9470 enfermos.
• Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce el primer día (0.338) con un valor de 9758 enfermos.

-Para [math]h=0.0001[/math]

Método de Euler

Método de Runge-Kutta

• Método de Euler: el momento en el que los enfermos son máximos se produce el primer día (0.3411) con un valor de 9465 enfermos.
• Método de Runge-Kutta: el momento en el que los enfermos son máximos se produce el primer día (0.3407) con un valor de 9493 enfermos.

Para este segundo caso, nos damos cuenta que al utilizar tamaño de paso [math]h=0.1[/math] , llegamos a un resultado sin sentido, pues nos da valores negativos y es imposible que una población tenga tales valores. En las demás gráficas, los resultados obtenidos no parecen cuadrar con la realidad. Esto podría ser debido a que la alta diferencia entre los valores iniciales [math](S _{0} , I_{0} )=(10000,40)[/math] no permiten un buen cálculo de su variación respecto del tiempo. El número inicial de infectados es muy elevado y al disminuir el tamaño de paso, el programa aproxima mejor la solución.

3 Resolución del problema completo suponiendo [math] a [/math] una función variable con el tiempo

Para el siguiente caso supondremos que el coeficiente [math] a [/math] (tasa de contagio en población sana) es una función dada por :

[math] a(t)= \frac{0.003}{(1+t)} [/math] Suponiendo además las condiciones iniciales [math] (S_{0},I_{0})=(1640,40) [/math] , vamos a dibujar las gráficas de infectados y sanos por el método de Heun, tomando un tamaño de paso [math] h=0.01 [/math]

Población infectada-Población sana para a=a(t)

%Heun. Coeficiente variable a(t)

%Datos
S0=input('Introducir valor inicial SO: ');
I0=input('Introducir valor inicial IO: ');
b=0.3;
c=0.01;
%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector solucion
y0=[S0;I0];
y=zeros(2,N+1);
y(:,1)=y0;
for i=1:N
     K1S=-(0.003/(1+t(i)))*y(1,i)*y(2,i);
     K2S=-(0.003/(1+t(i)+h))*(y(1,i)+K1S*h)*(y(2,i)+K1S*h);
     y(1,i+1)=y(1,i)+(h/2)*(K1S+K2S);
     K1I=(0.003/(1+t(i)))*y(1,i)*y(2,i)-b*y(2,i)-c*y(2,i);
     K2I=(0.003/(1+t(i)+h))*(y(1,i)+K1I*h)*(y(2,i)+K1I*h)-b*(y(2,i)+K1I*h)-c*(y(2,i)+K1I*h);
     y(2,i+1)=y(2,i)+(h/2)*(K1I+K2I);
end

%Tabla de resutados
[t',y(1,:)',y(2,:)']
%Graficos
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
title('Método de Heun');
legend('Poblacion sana','Población infectada','Location','best'); 
hold off

% Vector que contiene la variación de la población de infectados con el tiempo
I=y(2,:);
% Días de máximos infectados.
[fila,col]=find(I==max(max(I)));
% Valor máximo de infectados.
Diademaximo=(col-1)*h

Como se puede ver la población sana se infecta rápidamente y a partir del día 15 casi toda la población ha fallecido. Al segundo día [math](t=2,15)[/math] observamos que el número de infectados es de 994 individuos [math](I=994,3)[/math], siendo este el máximo que se alcanzará.

A medida que va pasando el tiempo el coeficiente [math] a(t) [/math] (tasa de contagio en población sana) va a ir decreciendo poco a poco hasta aproximarse a cero según podemos ver en el siguiente gráfico: [math] \lim_{t\to \infty } a(t)=0[/math]

a(t) para t=[0,40]

%tasa de contagio en poblacion sana a(t)

%Vector tiempo
h=input('Introducir tamaño de paso: ');
t0=0;
tN=40;
N=(tN-t0)/h;
t=t0:h:tN;

%Vector a
a=zeros(1,N+1);
for i=1:N+1
    a(i)=0.003/(1+t(i));
end

%Grafica
plot(t,a)

Poblaciones según la tasa de contagio a(t)

Ambas poblaciones parten del mismo inicio (para [math]t=0[/math] las poblaciones son iguales). Podemos observar que la población sana (susceptible de contraer la enfermedad) en el caso en el que la tasa de contagio en población sana [math](a)[/math] es constante [math](a=0.003)[/math] se ve reducida antes que en el caso en el que depende del tiempo [math](a=a(t))[/math].

[math]a(t)[/math] disminuye a la vez que transcurre el tiempo, por lo que en general las pendientes tanto de la población sana como de la infectada serán menores para este caso que para el caso en el que [math]a=0.003[/math]

4 Calibración del parámetro [math]a[/math] en base a una experiencia previa

Ahora vamos a intentar calibrar el párametro [math]a[/math] para que el número de infectados tenga el máximo en el quinto día, para ello escribimos el siguiente código:

%Aplicamos el Método de Heunn

%Datos
t0=0;
tN=40;
h=0.1;
t=t0:h:tN;
a=0.0005:0.0001:0.002;
b=0.3;
c=0.01;
N1=round((tN-t0)/h);
N2=length(a);
Tmax=zeros(1,N2);

%Vectores iniciales de Población e Infectados
y=zeros(2,N1+1);

%Condiciones iniciales
y(1,1)=1600;
y(2,1)=40;

%Bucle que nos calcula todas las posibles soluciones
for j=1:N2
    for i=1:N1
        K1S=-a(j)*y(1,i)*y(2,i);
        K2S=-a(j)*y(1,i)*y(2,i)+K1S*h;
        y(1,i+1)=y(1,i)+h/2*(K1S+K2S);
        
        K1I=y(2,i)*(a(j)*y(1,i)-(b+c));
        K2I=y(2,i)*(a(j)*y(1,i)-(b+c))+K1I*h;
        y(2,i+1)=y(2,i)+h/2*(K1I+K2I);
        
        if i>1 && y(2,i+1)<y(2,i) && y(2,i)>y(2,i-1)
            Tmax(j)=t(i);
        end
    end
end
distancia=(abs(Tmax-5));
[valormin,posmin]=min(distancia);
a=a(posmin)

Hemos encontrado la solución utilizando un bucle de tipo [math]for[/math] dos veces, esto nos permite obtener todas los valores de [math]I[/math], infectados, en el tiempo transcurrido entre [0,40] días para cualquier [math]a[/math] en el intervalo de [0.0005,0.002]. Una vez tenemos todas las posibilidades aplicamos una condición del tipo [math]if[/math] para localizar donde se sitúa el máximo de infectados en los distintos [math]a[/math] analizados, después vemos cuál de los máximos esta más próximo al quinto día, para ello generamos un vector [math]distancia[/math] donde hacemos el valor absoluto de la resta entre el valor obtenido y el quinto día tras esto conseguimos el valor mas cercano. Ahora analizamos el valor de [math]a[/math] para dicha posición obteniendo la solución pedida. En nuestro caso el valor pedida daba a = 8.0000e-04 que era el cuarto [math]a[/math] analizado en el bucle.