A mathematical model for the flu pandemic

A flu pandemic isn't only annoying: school children may get a few more free days, and traffic causes less air pollution. Furthermore, we can make a mathematical model of it, and experiment with it by changing the values of the parameters. It's always interesting to see which predictions are calculated by such a model.
I made for the Dutch mathematical magazine 'Onderzoek en Thuiscomputer' a simple model that can be extended in several ways. I also added a problem, which may be suitable as a (mini)puzzle for the readers of the magazine. Here it is.

Let:
t := time (in minutes); Δt: length of a small interval of time;
b0 := fraction of the infected in the beginning; this is a real number between 0 and 1 that indicates which fraction of the world population has been infected by the flu at the moment the governments are about to distribute on a large scale medicines against it;
s := duration (in minutes) of an individual infection (this infection ends either with recovery and immunity or (seldom) with death); as a first approximation, we suppose this duration is constant;
ΔT := length of the interval of time (in minutes) between two successive meetings of two persons in the world; again, as a first approximation, we suppose this very small interval length is constant;
ov := probability that during the meeting of an infected person with a person who is not yet infected or already immune the virus will hop from the infected person to the other one; again, as a first approximation, we suppose this probability is constant;
n := total number of persons in the world; again, as a first approximation, we suppose n is constant;
b(t) := fraction of the infected at point of time t;
g(t):=b(t-s); h(t):=b(t-2s); i(t):=b(t-3s); j(t):=b(t-4s); these are help functions I use in my pascal program to take into account, in the first five periods of s minutes since point of time t=0, that a person who is infected at a given point of time t, will be immune (or dead) s minutes later on;
c: =(2*ov)/(n*ΔT).

For t ≤ s we have (supposing immunity is negligible at the beginning of the pandemic):

b(t+ΔT) = b(t) + c*b(t)*(1-b(t))*ΔT

(, because, when two persons meet, the probability that one is infected and the other isn't, is equal to 2*b(t)*(1-b(t))).

So in approximation we have for t ≤ s and any small positive number Δt:

b(t+Δt) = b(t) + c*b(t)*(1-b(t))*Δt.

Problem: Which differential equation do you derive from this last equation? Solve the differential equation under the initial condition b(0)=b0. Prove that, if the differential equation would hold for t>s, b(t) would tend to 1 for t → ∞.

However, for t>s we have to take into account that a person who is infected at a given point of time t, will be immune (or dead) s minutes later on; we can do this using a pascal program like the following: (I tried to give the diverse parameters realistic values; they can varied anyway.)

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.

With the given values of the parameters, the output of the program is as follows (rounded to four places after the decimal point):

(At t=0 is b=0.01,) at t=6400 is b=0.0232, at t=12800 is b=0.0337, at t=19200 is b=0.0408, at t=25600 is b=0.0453, at t=32000 is b=0.0480.

So the fraction of the infected would in some three weeks increase from 1 percent to almost 5 percent, in a decreasing rate.

In order to refine the model, we may proceed as follows:
1) We partition the duration s of the infection into an incubation period s1 and a contagiousness (sickness) period s2.
2) We replace the supposition that the parameters s1, s2, ΔT and ov are constant with the supposition they are statistics with a normal distribution; here we approximate the normal distribution with mean m and standard deviation ss by a simpler triangle distribution from which we generate drawings with the following scheme:

For each meeting we draw a new value of the waiting time ΔT until the next meeting, a new value of ov, and, if the meeting leads to infection, new values of s1 and s2.
3) We distinguish between non-sick (number: n*(1-b)), incubating (number: inc) and immune (number: im). Only part of the non-sick are already immune, and another distinct part are still incubating.
4) To avoid capacity problems, we suppose the number of people in the world is only n=6000. We begin with 1 sick, 0 incubating and 0 immune. For ov we take as its mean value 0.1.
5) At the beginning of a meeting we add the new sick (and no more incubating) and the recovered (and henceforth immune) since the last meeting, and strike them away from the waiting line.
6) After the meeting we also draw a value from the homogeneous distribution on (0,1). If this value is smaller than 2*ov*b*(1-b-im/n-inc/n), we draw new values for s1 and s2 and place a sick in the waiting line at point of time t+s1, and a recovered at point of time t+s1+s2.

The new program is as follows:

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 see the number of momentaneous sick increases in only twenty-nine days to about 2500 (because of the high mean value 0.1 we chose for ov), and then decreases to 0 in eighteen days (because almost everybody is immune).
If we take 0.01 as the mean value of ov, the flu will not even spread. The first sick has already recovered before he could infect somebody.

Of course, teams at technological universities are already doing mathematical research into the course of pandemics more thoroughly.
But it's instructive to make simple models ourselves, experimenting with them and communicating about them. I'm curious to learn your comments.