[Wiskundige Modellering] Epidemiologische modellen

1182 days ago by Marnix.VanDaele

Deel II: Epidemiologische modellen

Hoofdstuk 2: Wiskundige modellering van besmettelijke ziekten

1. Het (herschaalde) SIR-model

\[
\left\{ \begin{array}{ll}
\displaystyle \frac{d\tilde{x}}{d\tau} = -R_0\, \tilde{x}\,\tilde{y} & \tilde{x}(0) = \frac{N-1}{N} \approx 1 \\[2mm]
\displaystyle\frac{d\tilde{y}}{d\tau} =R_0\, \tilde{x}\,\tilde{y} -\tilde{y} \qquad & \tilde{y}(0)=\frac{1}{N} \approx 0 \\[2mm]
\displaystyle\frac{d\tilde{z}}{d\tau} = \tilde{y} & \tilde{z}(0)=0 \,.
\end{array} \right.\]

SIR-model:  Verloop van de oplossing in functie van de tijd voor een specifieke waarde van $R_0$.

# SIR - verloop van de oplossing in functie van de tijd voor een specifieke waarde van R0 # startwaarde: (0.999,0.001,0) # reset() from sage.calculus.desolvers import desolve_odeint from sage.calculus.functions import jacobian x,y,z=var('x,y,z') R0=2; f=[-R0*x*y,R0*x*y-y,y] sr=srange(0,20,0.1); sol=desolve_odeint(f,[0.999,0.001,0],sr,[x,y,z]) px=line(zip(sr,sol[:,0]),color='black',legend_label=r'$\tilde{x}(\tau)$'); py=line(zip(sr,sol[:,1]),color='red',legend_label=r'$\tilde{y}(\tau)$'); pz=line(zip(sr,sol[:,2]),color='blue',legend_label=r'$\tilde{z}(\tau)$'); show(px+py+pz,axes_labels=[r'$\tau$',''],aspect_ratio=20, title='SIR-model met $R_0=3$') 
       

Het (herschaalde) SIR-model: resterende fractie van vatbare personen als functie van $R_0$.

# SIR - Het (herschaalde) SIR-model: resterende fractie van vatbare personen als functie van R0 # var('R0') implicit_plot(R0+log(x)/(1-x),(R0,1,5),(x,0,1),axes_labels=[r'$R_0$',r'$\tilde{x}_\infty$'],aspect_ratio=4,color='blue') 
       

Het herschaalde SIR-model:  plot van de fractie personen die de ziekte heeft opgelopen na afloop van epidemie als functie van $R_0$

# SIR - plot van fractie personen die ziekte heeft opgelopen na afloop van epidemie als functie van R0 # implicit_plot(R0+log(1-x)/(x),(R0,1,5),(x,0,1),axes_labels=[r'$R_0$',r'$\tilde{z}_\infty$'],aspect_ratio=4,color='blue') 
       

Het (herschaalde) SIR-model: Plot van de maximale graad van infectie en corresponderde maximale graad van vatbaarheid als functie van $R_0$

# SIR - plot van maximale graad van infectie en corresponderde maximale graad van vatbaarheid als functie van R0 # # p1=plot(1-(1+log(R0))/R0,(R0,1,5),color='red',legend_label='$1-(1+\log(R_0))/R_0$') p2=plot(1/R0,(R0,1,5),color='black',axes_labels=['$R_0$',''],legend_label='$1/R_0$') show(p1+p2) 
       

Het herschaalde SIR-model: verloop van de oplossing in het ($\tilde{x}$, $\tilde{y}$)-vlak voor verschillende waarden van $R_0$. 

# SIR - verloop van de oplossing in het (\tilde{x},\tilde{y})-vlak voor # verschillende waarden van R0 # startwaarde: (0.999,0.001,0) # # reset() from sage.calculus.desolvers import desolve_odeint from sage.calculus.functions import jacobian x,y,z=var('x,y,z') R0=8; sr=srange(0,50,0.05); g = Graphics() g=g+ line([[1, 0],[0,1]],color='black') while R0>0.9: f=[-R0*x*y,R0*x*y-y,y] sol=desolve_odeint(f,[0.999,0.001,0],sr,[x,y,z]); p=line(zip(sol[:,0],sol[:,1]),color=hue(R0/5+0.5),legend_label=R0); g=g+p R0=R0-1 g=g+ parametric_plot([1/x,1-(1+log(x))/x],(x,1,10),color='black',linestyle='dashed'); g.show(dpi=100,axes_labels=[r'$\tilde{x}$',r'$\tilde{y}$'],aspect_ratio=0.75) 
       

2. Het herschaalde SIRS-model

\[
\left\{ \begin{array}{ll}
\displaystyle \frac{d\tilde{x}}{d\tau} = -\frac{\beta}{\gamma}\, \tilde{x}\,\tilde{y} + \frac{\omega}{\gamma}\,\tilde{z} \qquad & \tilde{x}(0) = \frac{N-1}{N} \approx 1 \\[2mm]
\displaystyle \frac{d\tilde{y}}{d\tau} =\frac{\beta}{\gamma}\, \tilde{x}\,\tilde{y} -\tilde{y} & \tilde{y}(0)=\frac{1}{N} \approx 0 \\[2mm]
\displaystyle \frac{d\tilde{z}}{d\tau} = \tilde{y} - \frac{\omega}{\gamma} \tilde{z}& \tilde{z}(0)=0 \,.
\end{array} \right. \]


Het (herschaalde) SIRS-model: studie van de stabiliteit van de evenwichtspunten

# SIRS - studie van de stabiliteit van de evenwichtspunten reset; from sage.calculus.functions import jacobian x,y,z,R0,om,omega,gamma=var('x,y,z,R0,om,omega,gamma') # om=omega/gamma; f=[-R0*x*y+om*z,R0*x*y-y,y-om*z] g=[f[0],f[1],x+y+z-1]; # opleggen dan x+y+z=1 sol=solve(g,{x,y,z}) J=jacobian([f[0].subs(z=1-x-y),f[1].subs(z=1-x-y)],[x,y]) pretty_print('jacobiaan:',J) pretty_print('evenwichtspunten:',sol[0],'en ', sol[1]) J0=J.subs(sol[0]) print('eerste evenwichtspunt: eigenwaarden van jacobiaan') pretty_print(J0.eigenvalues()) J1=simplify(J.subs(sol[1])) cp=J1.characteristic_polynomial() print('tweede evenwichtspunt: karakteristieke veelterm van jacobiaan') pretty_print(simplify(cp(x))) coef=cp(x).coefficients(x) pretty_print('som eigenwaarden=',-factor(coef[1][0]/coef[2][0])); pretty_print('product eigenwaarden=',factor(coef[0][0]/coef[2][0])); 
       
eerste evenwichtspunt: eigenwaarden van jacobiaan
tweede evenwichtspunt: karakteristieke veelterm van jacobiaan

                                
                            
eerste evenwichtspunt: eigenwaarden van jacobiaan
tweede evenwichtspunt: karakteristieke veelterm van jacobiaan

                                

Het (herschaalde) SIRS-model: Verloop van de oplossing in functie van de tijd voor specifieke waarden van $R_0=\beta/\gamma$ en $om= \omega/\gamma$

# SIRS - verloop van de oplossing in functie van de tijd voor een # specifieke waarde van R0 en om # startwaarde: (0.999,0.001,0) # from sage.calculus.desolvers import desolve_odeint from sage.calculus.functions import jacobian x,y,z=var('x,y,z') R0=3; om=0.20; # om= omega/gamma f=[-R0*x*y+om*z,R0*x*y-y,y-om*z] sr=srange(0,20,0.1); sol=desolve_odeint(f,[0.999,0.001,0.0],sr,[x,y,z]) px=line(zip(sr,sol[:,0]),color='black',legend_label=r'$\tilde{x}(\tau)$'); py=line(zip(sr,sol[:,1]),color='red',legend_label=r'$\tilde{y}(\tau)$'); pz=line(zip(sr,sol[:,2]),color='blue',legend_label=r'$\tilde{z}(\tau)$'); show(px+py+pz,axes_labels=[r'$\tau$',''],title=r'SIRS-model met $R_0=\beta/\gamma=3$ en $\omega/\gamma=0.20$ met beginwaarden $(0.999, 0.001, 0)$') 
       

Het (herschaalde) SIRS-model: Verloop van de oplossing in het ($\tilde{x}$, $\tilde{y}$)-vlak voor verschillende waarden van $R_0=\beta/\gamma$ en vaste $\omega/\gamma$

# SIRS - verloop van de oplossing in het (\tilde{x}, \tilde{y})-vlak # voor verschillende waarden van R0=beta/gamma en vaste omega/delta # startwaarde: (0.999,0.001,0) of (0.3,0.5, 0.2) # # #reset() from sage.calculus.desolvers import desolve_odeint from sage.calculus.functions import jacobian x,y,z=var('x,y,z') sp=[0.999,0.001,0] # startpunt R0=20; om=0.2 # om = omega/gamma sr=srange(0,50,0.05); g = Graphics() gg = Graphics() l_ymax=[] l_ymaxR0=[] g=g+ line([[1, 0],[0,1]],color='black') while R0>0.9: f=[-R0*x*y+om*z,R0*x*y-y,y-om*z] sol=desolve_odeint(f,sp,sr,[x,y,z]); if (abs(R0-round(R0))<0.1) & (R0<9): p=line(zip(sol[:,0],sol[:,1]),color=hue(R0/5+0.5),legend_label=R0); g=g+p ind=0 s=sol[:,1] leng=len(s) while (s[ind]>s[ind+1]) & (ind<leng-2): ind=ind+1; ss=s[ind:] if ind==0: l_ymax.append([1/R0,max(sol[:,1])]) if R0<=8: l_ymaxR0.append([R0,max(sol[:,1])]) else: if R0<=8: l_ymaxR0.append([R0,sp[1]]) if (abs(R0-round(R0))<0.1) & (R0<9) : g=g+ line([[1/R0, 0],[1/R0, max([max(ss),(1-1/R0)*om/(1+om)])]],color=hue(R0/5+0.5),linestyle='dotted') # verticale stippellijn R0=R0-0.2 # R0=8 # l=[] # lijn van minimale x # while R0>0.9: # f=[-R0*x*y+om*z,R0*x*y-y,y-om*z] # sol=desolve_odeint(f,sp,sr,[x,y,z]); # xm=min(sol[:,0]) # l.append([xm,om*(1-xm)/(R0*xm+om)]) # R0=R0-0.25 # g=g+line(l,color='gray',linestyle='dashed') # minimale x g=g+line(l_ymax,color='gray',linestyle='dashed') # maximale y g=g+ parametric_plot([1/x,(1-1/x)*om/(1+om)],(x,1,8),color='black',linestyle='dotted'); #gg=gg+ parametric_plot([x,(1-1/x)*om/(1+om)],(x,1,8),linestyle='dotted',legend_label=r'$\tilde{y}_\infty$',legend_label=om); # evenwichtspunten gg=gg+ line(l_ymaxR0,color='gray',linestyle='dashed',legend_label=r'$\tilde{y}_{\rm max}$') # maximale y g.show(dpi=100,axes_labels=[r'$\tilde{x}$',r'$\tilde{y}$'],aspect_ratio=0.75) gg.show(dpi=100,axes_labels=[r'$R_0$',r'${\ }$'],aspect_ratio=10) 
       

Het (herschaalde)  SIRS-model: Maximale waarde $\tilde{y}_{\rm max}$ als functie van $R_0=\beta/\gamma$ voor verschillende waarden van $\omega/\delta$

# SIRS - maximale waarde \tilde{y}_max als functie van R0=beta/gamma # voor verschillende waarden van omega/delta # startwaarde: (0.999,0.001,0) of (0.3,0.5, 0.2) # reset() from sage.calculus.desolvers import desolve_odeint from sage.calculus.functions import jacobian x,y,z=var('x,y,z') sp=[0.999,0.001,0] # startpunt om=0.0 # om = omega/gamma sr=srange(0,50,0.05); gg = Graphics() ggg = Graphics() while om<5: R0=20 l_ymax=[] l_ymaxR0=[] while R0>0.9: f=[-R0*x*y+om*z,R0*x*y-y,y-om *z] sol=desolve_odeint(f,sp,sr,[x,y,z]); ind=0 s=sol[:,1] leng=len(s) while (s[ind]>s[ind+1]) & (ind<leng-2): ind=ind+1; ss=s[ind:] if ind==0: l_ymax.append([1/R0,max(sol[:,1])]) if R0<=8: l_ymaxR0.append([R0,max(sol[:,1])]) else: if R0<=8: l_ymaxR0.append([R0,sp[1]]) R0=R0-0.2 gg=gg+ line(l_ymaxR0,linestyle='dashed',legend_label=om,color=hue(om/5+0.5)) # maximale y ggg=ggg+ parametric_plot([x,(1-1/x)*om/(1+om)], (x,1,8),legend_label=om,linestyle='dotted',color=hue(om/5+0.5)); # # evenwichtspunten om=om+0.5 gg.show(dpi=100,axes_labels=[r'$R_0$',r'$\tilde{y}_{\rm max}$'],aspect_ratio=5) ggg.show(dpi=100,axes_labels=[r'$R_0$',r'$\tilde{y}_{\infty}$'],aspect_ratio=8) 
       

3. Het (herschaalde) SIR-model met demografie

 

\[
\left\{ \begin{array}{l}
\displaystyle \frac{d\tilde{x}}{d\tau} = \frac{\mu}{\gamma} -\frac{\beta}{\gamma}\, \tilde{x}\,\tilde{y} - \frac{\mu}{\gamma}\,\tilde{x} \\[2mm]
\displaystyle \frac{d\tilde{y}}{d\tau} =\frac{\beta}{\gamma}\, \tilde{x}\,\tilde{y} - \left(1+\frac{\mu}{\gamma}\right)\tilde{y} \\[2mm]
\displaystyle \frac{d\tilde{z}}{d\tau} = \tilde{y} - \frac{\mu}{\gamma} \tilde{z} \,.
\end{array} \right.\]



Het (herschaalde) SIR-model met demografie - verloop van de oplossing in functie van de tijd voor een specifieke waarde van beta/gamma en mu/gamma

# SIR met demografie - verloop van de oplossing in functie van de tijd voor # een specifieke waarde van beta/gamma en mu/gamma # startwaarde: (0.999,0.001,0) # # reset() from sage.calculus.desolvers import desolve_odeint from sage.calculus.functions import jacobian x,y,z=var('x,y,z') # A= mu/gamma # B= beta/gamma A=0.1 B=5 # R0=beta/(gamma+mu)=1/(gamma/beta+(mu/gamma)*(gamma/beta))=1/(1/B+A/B))=B/(1+A) R0=B/(1+A) pretty_print(R0) f=[-B*x*y+A*(1-x),B*x*y-(1+A)*y,y-A*z] sr=srange(0,20,0.1); sol=desolve_odeint(f,[0.999,0.001,0],sr,[x,y,z]) px=line(zip(sr,sol[:,0]),color='black',legend_label=r'$\tilde{x}(\tau)$'); py=line(zip(sr,sol[:,1]),color='red',legend_label=r'$\tilde{y}(\tau)$'); pz=line(zip(sr,sol[:,2]),color='blue',legend_label=r'$\tilde{z}(\tau)$'); show(px+py+pz,axes_labels=[r'$\tau$',''],title=r'SIR-model met demografie met $\beta/\gamma=0.1$ en $\mu/\gamma=5$') 
       

                                
                            

                                

Het (herschaalde) SIR-model met demografie: Verloop van de oplossing in het  ($\tilde{x}$, $\tilde{y}$)-vlak voor verschillende waarden van $\beta/\gamma$ en $\mu/\gamma= 0.2$

S# SIR met demografie - verloop van de oplossing in het ($\tilde{x}$, $\tilde{y}$)-vlak voor verschillende waarden van \beta/gamma en \mu/gamma= 0.2 vast # startwaarde: (0.999,0.001,0) of (0.3,0.5,0.2) # A=0.2 = mu/gamma # # A=0.2; B=8; sr=srange(0,60,0.05); g = Graphics() g=g+ line([[1, 0],[0,1]],color='black') while B>0.2: f=[-B*x*y+A*(1-x),B*x*y-(1+A)*y,y-A*z] sol=desolve_odeint(f,[0.3,0.5,0.2],sr,[x,y,z]); p=line(zip(sol[:,0],sol[:,1]),color=hue(B/5+0.5)); g=g+p B=B-1 g=g+ parametric_plot([1/x,(1-1/x)*A/(1+A)],(x,1,10),color='black',linestyle='dotted'); g.show(dpi=100,axes_labels=[r'$\tilde{x}$',r'$\tilde{y}$']) 
       

Het (herschaalde) SIR-model met demografie: studie van de stabiliteit van evenwichtspunten

# studie van stabiliteit van evenwichtspunten in SIR model met demografie reset(); x,y,z,mu,beta,gamma,R0=var('x,y,z,mu,beta,gamma,R0') A= mu/gamma B= beta/gamma f=[-B*x*y+A*(1-x),B*x*y-(1+A)*y,y-A*z] sol=solve(f,[x,y,z]) #pretty_print(sol) ziektevrij=sol[0] pretty_print('ziektevrij evenwichtspunt:', ziektevrij) g=[-B*x*y+A*(1-x),B*x*y-(1+A)*y] jac=jacobian(g,[x,y]) J1=jac.subs(ziektevrij) tr=factor(J1.trace()) dt=factor(J1.det()) pretty_print('som van eigenwaarden = spoor = ',tr) pretty_print('product van eigenwaarden = det =',dt) endem=sol[1] pretty_print('endemisch evenwichtspunt:', endem) g=[-B*x*y+A*(1-x),B*x*y-(1+A)*y] jac=jacobian(g,[x,y]) J1=jac.subs(endem) tr=factor(J1.trace()) dt=factor(J1.det()) pretty_print('som van eigenwaarden = spoor = ',tr) pretty_print('product van eigenwaarden = det =',dt) #discr=factor(tr^2-4*dt) #pretty_print('discr=',discr.factor())