Een grieppandemie heeft ook leuke kanten: de scholieren krijgen wellicht een paar dagen vrij, en er is minder fijnstof door autoverkeer. En, nog mooier, je kunt
er een wiskundig model voor maken, en daarmee experimenteren door de waarden van de parameters te variëren. Het is altijd spannend om te bekijken welke
voorspellingen zo'n model berekent.
Ik heb speciaal voor Onderzoek & Thuiscomputer een eenvoudig model gemaakt dat op verschillende manieren verfijnd kan worden. Er is ook een opgave aan
verbonden die misschien geschikt is als (mini)puzzel voor dit interessante wiskundepuzzelblad. Hier komt ie.
Stel:
t := tijd (in minuten); Δt: lengte van een zeer klein tijdsinterval;
b0 := besmettingsgraad 'aan het begin'; dit is een getal tussen 0 en 1 dat aangeeft welk deel van de wereldbevolking de griep onder de leden heeft op
het moment dat men besluit om op grote schaal griepwerende middelen te verspreiden;
s := duur (in minuten) van een individuele besmetting (deze eindigt ofwel met genezing en immuniteit ofwel (zelden) met de dood); deze duur wordt hier, als een
eerste benadering, constant verondersteld;
ΔT := tijdsinterval (in minuten) tussen twee opeenvolgende ontmoetingen van twee wereldburgers; deze zeer korte tijd wordt hier eveneens, als een eerste
benadering, constant verondersteld;
ov := kans op virusoverdracht bij ontmoeting van een besmet individu met een onbesmet (of reeds immuun) individu; ook deze kans wordt hier, als een eerste
benadering, constant verondersteld;
n := aantal wereldburgers; ook n wordt in eerste benadering constant verondersteld;
b(t) := besmettingsgraad op tijdstip t;
g(t):=b(t-s); h(t):=b(t-2s); i(t):=b(t-3s); j(t):=b(t-4s); dit zijn hulpfuncties om voor de eerste vijf perioden van s minuten sinds tijdstip 0 in mijn
pascalprogramma te verdisconteren dat een persoon die op een gegeven moment besmet wordt, s minuten later immuun is (of overleden);
c: =(2*ov)/(n*ΔT).
Voor t ≤ s geldt (uitgaande van de veronderstelling dat immuniteit aan het begin van de pandemie nog geen rol speelt):
b(t+ΔT) = b(t) + c*b(t)*(1-b(t))*ΔT
(,want bij een ontmoeting tussen twee wereldburgers is de kans dat er één besmet is en de ander niet, gelijk is aan 2*b(t)*(1-b(t))).
Dus geldt bij benadering voor t ≤ s en elk klein positief getal Δt:
b(t+Δt) = b(t) + c*b(t)*(1-b(t))*Δt.
Opgave: Stel de differentiaalvergelijking op die hoort bij deze differentievergelijking; los hem op onder beginvoorwaarde b(0)=b0; bewijs dat b(t) naar 1 zou naderen voor t → ∞ als deze differentiaalvergelijking ook voor t>s zou gelden.
Voor t>s moeten we echter verdisconteren dat een persoon die op een gegeven moment besmet wordt, s minuten later immuun is (of overleden); dit kunnen we doen met een pascalprogramma als het volgende: (waarin ik de diverse parameters realistische waarden heb proberen te geven; deze kunnen echter gevarieerd worden.)
program pandemie;
var s,t:integer; n,b0,DT,ov,c,b,g,h,i,j:real;
begin
b0:=0.01;s:=6400;DT:=1/40000000;ov:=0.01;n:=6000000000;
c:=(2*ov)/(n*DT);
b:=b0;
for t:=1 to s do
begin
b:=b+c*b*(1-b)
end;
writeln(t:7,b:13:10);readln;
g:=b0;
for t:=s+1 to 2*s do
begin
g:=g+c*g*(1-g);
b:=b+c*(b*(1-b)-g*(1-g))
end;
writeln(t:7,b:13:10);readln;
h:=b0;
for t:=2*s+1 to 3*s do
begin
h:=h+c*h*(1-h);
g:=g+c*(g*(1-g)-h*(1-h));
b:=b+c*(b*(1-b)-g*(1-g))
end;
writeln(t:7,b:13:10);readln;
i:=b0;
for t:=3*s+1 to 4*s do
begin
i:=i+c*i*(1-i);
h:=h+c*(h*(1-h)-i*(1-i));
g:=g+c*(g*(1-g)-h*(1-h));
b:=b+c*(b*(1-b)-g*(1-g))
end;
writeln(t:7,b:13:10);readln;
j:=b0;
for t:=4*s+1 to 5*s do
begin
j:=j+c*j*(1-j);
i:=i+c*(i*(1-i)-j*(1-j));
h:=h+c*(h*(1-h)-i*(1-i));
g:=g+c*(g*(1-g)-h*(1-h));
b:=b+c*(b*(1-b)-g*(1-g))
end;
writeln(t:7,b:13:10);readln;
end.
Bij de gegeven waarden van de parameters is de output van dit programma als volgt (afgerond op vier decimalen na de decimaalpunt):
(Op t=0 is b=0.01,) op t=6400 is b=0.0232, op t=12800 is b=0.0337, op t=19200 is b=0.0408, op t=25600 is b=0.0453, op t=32000 is b=0.0480.
Dus dan zou de besmettingsgraad in ruim drie weken tijd in een afnemend tempo oplopen van 1 percent tot bijna 5 percent.
Om het model te verfijnen, kunnen als volgt te werk gaan:
1) We verdelen de besmettingsduur s in een incubatietijd s1 en een besmettelijkheidstijd (ziekzijntijd) s2.
2) We vervangen de veronderstelling dat de parameters s1, s2, ΔT en ov constant zijn door de aanname dat het normaal
verdeelde stochasten zijn; hierbij benaderen we de normale verdeling met gemiddelde m en standaardafwijking ss door een driehoeksverdeling
waaruit de trekkingen eenvoudig worden gegenereerd met het volgende schema:
Bij elke ontmoeting doen we nieuwe trekkingen voor de wachttijd ΔT tot de volgende ontmoeting, voor ov, en, als de ontmoeting een infectie tot
gevolg heeft, voor s1 en s2.
3) We maken onderscheid tussen niet-zieken (aantal: n*(1-b)), incuberenden (aantal: inc) en immunen (aantal: im). Slechts een deel van de niet-zieken is al
immuun, en een ander deel is nog in de incubatietijd.
4) Om capaciteitsproblemen te vermijden, stellen we het aantal wereldburgers op n=6000. We beginnen met 1 zieke, 0 incuberenden en 0 immunen.
Voor ov nemen we een gemiddelde waarde van 0.1.
5) Bij aanvang van een ontmoeting verwerken we de ziekgewordenen (en dus niet meer incuberenden) en de genezenen (en dus nu immunen)
sinds de vorige ontmoeting, en schrappen ze uit de wachtrij.
6) Na elke ontmoeting doen we een trekking uit de homogene verdeling op (0,1). Als die kleiner is dan 2*ov*b*(1-b-im/n-inc/n), doen we een trekking
voor s1 en een voor s2 en zetten een nieuwe zieke in de wachtrij op tijdstip t+s1, en een genezene op tijdstip t+s1+s2.
Het nieuwe programma ziet er als volgt uit:
program pandemic;
type rij = array[1..maxint,1..2] of real;
var r,t,s1,s2,dt0,dt,ov,b,z,im,inc:real; k,l,n:integer; p:rij;
function trekking(m,ss:real):real;
var r,a:real;
begin
r:=random;
if not r≥0.5 then a:=sqrt(2*r)-1 else a:=1-sqrt(2-2*r);
trekking:=m+a*ss
end {trekking};
procedure initialiseer(var p:rij);
var i:integer;
begin
for i:=1 to maxint do
begin p[i,1]:=1000000; p[i,2]:=0 end
end {initialiseer};
procedure voegtoe(k:integer; t:real; d:integer; var p:rij);
var i,plaats:integer; klaar:boolean;
begin
i:=k; klaar:=false;
repeat
if i=1 then
begin klaar:=true; if not p[1,1]≥t then plaats:=2 else plaats:=1 end
else if not p[i-1,1]≥t then begin plaats:=i; klaar:=true end else i:=i-1
until klaar;
for i:=k downto plaats+1 do begin p[i,1]:=p[i-1,1]; p[i,2]:=p[i-1,2] end;
p[plaats,1]:=t; p[plaats,2]:=d
end {voegtoe};
procedure verwerk(var p:rij; t:real; var l:integer; var z,im,inc:real);
var i:integer; klaar:boolean;
begin
l:=0; klaar:=false; i:=1;
while not klaar do
begin
if not p[i,1]≥t then begin z:=z+p[i,2]; l:=l+1; if not p[i,2]≥0 then im:=im+1 else inc:=inc-1 end else klaar:=true;
i:=i+1
end
end {verwerk};
procedure vergeet(l:integer; var k:integer; var p:rij);
var i:integer;
begin
for i:=1 to k-l do
begin p[i,1]:=p[i+l,1]; p[i,2]:=p[i+l,2] end;
for i:=k-l+1 to k do
begin p[i,1]:=10000000; p[i,2]:=0 end;
k:=k-l
end {vergeet};
begin
t:=0;
n:=6000;
initialiseer(p);
z:=1; im:=0; inc:=0; b:=z/n;
k:=1; s2:=trekking(3000,1500);
voegtoe(k,s2,-1,p);
dt0:=1/(((10*n)/2)/(24*60));
repeat
dt:=trekking(dt0,dt0/2);
t:=t+dt;
verwerk(p,t,l,z,im,inc); b:=z/n;
vergeet(l,k,p);
writeln(t:13:10,' *',round(z):7,' *',round(inc):7,' *',round(im):7);
ov:=trekking(0.1,0.05);
r:=random;
if not r≥2*ov*b*(1-b-im/n-inc/n) then
begin
inc:=inc+1;
s1:=trekking(3000,1500); s2:=trekking(6000,3000);
k:=k+1; voegtoe(k,t+s1,1,p);
k:=k+1; voegtoe(k,t+s1+s2,-1,p)
end;
until false
end.
We zien dat het aantal gelijktijdig zieken bij deze hoge gemiddelde waarde van ov (0.1) in slechts negenentwintig dagen oploopt tot ongeveer 2500, en daarna in
achttien dagen afneemt tot 0 (omdat bijna iedereen immuun is).
Voor ov gemiddeld 0.01 komt de griep niet eens van de grond. De eerste zieke is al weer beter voordat hij iemand besmet heeft.
Uiteraard wordt wiskundig onderzoek naar het verloop van pandemieën op de technische universiteiten al lang en grondiger
en in teamverband uitgevoerd.
Het is niettemin leerzaam als we zelf ook eenvoudige modellen maken, ermee experimenteren, en erover communiceren.
Ik ben benieuwd naar uw commentaar.