title: Een korte geschiedenis van wiskundige populatiedynamica. , 164 p., 2022, 979-10-343-6356-8. &#x27E8;hal-03551687&#x27E9;
text: George Udny Yule werd in 1871 in Schotland geboren. Zijn vader bekleedde een hoge functie in het Britse bestuur in India. Op 16-jarige leeftijd begon Yule te studeren aan het University College in Londen om ingenieur te worden. In 1892 veranderde hij van richting en deed hij een jaar onderzoek in Bonn onder supervisie van de natuurkundige Heinrich Hertz, die een paar jaar eerder het bestaan van elektromagnetische golven had aangetoond. Toen Yule terugkeerde naar Engeland, bood Karl Pearson hem een positie aan als universitair docent in de toegepaste wiskunde aan het University College. Yule begon zich, in navolging van Pearson, te richten op statistiek. In 1911 publiceerde hij Een inleiding tot de theorie van de statistiek, dat 14 maal herdrukt werd. Het jaar daarop verhuisde hij naar de universiteit van Cambridge. Zijn onderzoek had betrekking op theoretische aspecten van de statistiek, maar ook op toepassingen in de landbouw en de epidemiologie. In 1922 werd hij lid van de Royal Society. In 1924 publiceerde Yule een artikel getiteld Een wiskundige evolutietheorie gebaseerd op de conclusies van Dr. J. C. Willis. Willis was een collega van de Royal Society die in 1922 een boek had gepubliceerd met de titel Leeftijd en gebied, een studie in geografische distributie en oorsprong van soorten. Hij had de verdeling van soorten over verschillende geslachten in de classificatie van planten en dieren bestudeerd. Uit de gegevens die hij had verzameld, bleek dat de meeste geslachten slechts één soort bevatten, dat steeds minder geslachten een groter aantal soorten bevatten en dat er nog een paar geslachten waren met een groot aantal soorten. In tabel 15.1 staan de 84 Figuur 15.1: Yule (1871-1951) gegevens betreffende slangen, hagedissen en twee keverfamilies (de Chrysomelidae en de Cerambycinae). De 1.580 toen bekende soorten hagedissen waren ingedeeld in 259 geslachten, waarvan 105 geslachten slechts één soort bevatten, 44 slechts twee soorten, 23 slechts drie soorten, enz. en twee geslachten met meer dan honderd soorten. Voor andere families van dieren en planten had de verdeling van de geslachten naar het aantal soorten dat zij bevatten een zeer vergelijkbare vorm. Yule stelde voor dat Willis zou proberen zijn gegevens uit te zetten in een grafiek met logaritmische schalen. Dit leverde een opmerkelijk resultaat op (figuur 15.2): de logaritme van het aantal Qn geslachten die n soorten bevatten, neemt min of meer lineair af met log(n). Met andere woorden, er zijn constanten α > 0 en β > 0 zodanig dat Qn ≈ α n−β : de verdeling volgt een „machtswet”. In zijn artikel van 1924 zocht Yule naar een wiskundig evolutiemodel dat zo’n statistische verdeling kon verklaren. Figuur 15.2: Het aantal geslachten als functie van het aantal soorten dat ze bevatten, met decimale logaritmische schalen. Gegevens voor de Chrysomelidae. Om de schommelingen bij een groot aantal soorten n af te vlakken, werden de geslachten geteld voor een interval van n-waarden, zoals in tabel 15.1. Het gemiddelde aantal geslachten voor één enkele waarde van n kan dus minder dan 1 zijn.
  
  aantal genus 2 10 1 10 0 10 aantal soorten −1 10 0 10 1 10 2 10 Hoofdstuk 15 85 Tabel 15.1: Gegevens verzameld door Willis. Aantal soorten Chrysomelidae Aantal geslachten Cerambycinae Slangen 1 2 3 4 5 6 7 8 9 10 11-20 21-30 31-40 41-50 51-60 61-70 71-80 81-90 91-100 101- 215 90 38 35 21 16 15 14 5 15 58 32 13 14 5 8 7 7 3 16 469 152 82 61 33 36 18 17 14 11 74 21 15 8 4 3 0 1 1 4 131 35 28 17 16 9 8 8 9 4 10 12 3 1 0 0 1 0 1 0 105 44 23 14 12 7 6 4 5 5 17 9 3 2 0 1 0 0 0 2 totaal 627 1024 293 259 Hagedissen
  
  Daartoe bedacht hij eerst een stochastisch model1 met continue tijd voor de groei van het aantal soorten binnen één geslacht (figuur 15.3a). Uitgaande van slechts één soort op het tijdstip t = 0, veronderstelde hij dat de kans voor een soort om door mutatie het leven te schenken aan een nieuwe soort van hetzelfde geslacht gedurende een „klein” tijdsinterval dt (op de tijdschaal van de evolutie) gelijk was aan r dt met r > 0. Zij pn (t) de waarschijnlijkheid dat er n soorten zijn op tijdstip t (n is een geheel getal maar t is een reëel getal). Om pn (t + dt) te berekenen, beschouwde Yule verschillende gevallen: • als er op het tijdstip t n − 1 soorten zijn, heeft elke soort een kans r dt om tussen t en t + dt één nieuwe soort voort te brengen; in de limiet 1 McKendrick (zie hoofdstuk 16) was reeds begonnen met het bestuderen van dergelijke modellen in de populatiedynamica in een in 1914 gepubliceerd artikel.
  
  86 1 1 1 3 n=1 3 1 4 n=2 n=3 2 t 0 n=4 t
  
  Figuur 15.3: (a) Een simulatie van de evolutie van het aantal soorten binnen één geslacht. Soort 1 genereert soorten 2 en 3. Soort 3 genereert soort 4. (b) De kans pn (t) dat er n soorten van hetzelfde geslacht zijn op tijdstip t, voor 1 n 4. dt → 0 zullen er n soorten zijn op het tijdstip t + dt met een kans (n − 1) r dt; • als er n soorten zijn op tijdstip t, zullen er n + 1 soorten zijn op tijdstip t + dt met een kans van n r dt. Dus pn (t) wordt gegeven door het volgende stelsel differentiaalvergelijkingen d p1 = −r p1, dt d pn = (n − 1) r pn−1 − n r pn dt (15.1) (15.2) voor alle n 2. Uit de eerste vergelijking volgt p1 (t) = e−r t omdat p1 (0) = 1. Men kan aantonen dat de oplossing van de tweede vergelijking met pn (0) = 0 gegeven wordt door pn (t) = e−r t (1 − e−r t )n−1 (15.3) voor alle n 2 (figuur 15.3b). Dus op een vast tijdstip t is de verdeling van de kansen (pn (t))n1 „geometrisch” en is de verhouding tussen twee opeenvolgende termen gelijk aan 1 − e−r t. We merken immers eerst op dat vergelijking (15.2) equivalent is met d pn en r t = (n − 1) r pn−1 en r t, dt (15.4) waaruit we achtereenvolgens p2 (t), p3 (t).... kunnen berekenen. We krij- p2 (t) = e−r t (1 − e−r t ), dan p3 (t) = e−r t (1 − e−r t )2, wat formule (15.3) suggereert voor de algemene oplossing. Men kan tenslotte nagaan dat deze formule een oplossing is van vergelijking (15.4). Yule leidde ook uit formule (15.3) af dat het verwachte aantal soorten +∞ exponentieel toeneemt met de tijd: ∑ n pn
  
  (
  
  t) =
  
  er 
  
  t
  
  . n=1
  
  We merken immers eerst op dat
  
  voor |x| < 1, +∞ ∑ nx n−1 n=1 Dan d +∞ n d = ∑ x = dx dx n=0 +∞ +∞ n=1 n=1 1 1−x = 1. (1 − x)2 ∑ n pn (t) = e−r t ∑ n(1 − e−r t )n−1 = er t. In het bijzonder, als T de verdubbelingstijd is, gedefinieerd door erT = 2, dan is de kansverdeling (pn (t))n1 van het aantal soorten op het tijdstip t = T 1 meetkundig met een verhouding 1/2: 21, 14, 18, 16... Op het tijdstip t = kT is k het meetkundig met een verhouding 1 − 1/2 en p1 (kT ) = 1/2k. Yule beschouwde vervolgens, parallel met de groei van het aantal soorten dat tot hetzelfde geslacht behoort, een soortgelijk proces als gevolg van grotere mutaties die leiden tot het ontstaan van nieuwe geslachten. Zij s dt de kans voor een bestaand geslacht om een nieuw geslacht te genereren gedurende een klein tijdsinterval dt. Als we ervan uitgaan dat er op het tijdstip t = 0 slechts één geslacht is, wordt het verwachte aantal geslachten op het tijdstip t gegeven door est. Het gemiddeld aantal geslachten dat per tijdseenheid op tijdstip t ontstaat is de afgeleide s est. In de limiet2 t → +∞, is het gemiddelde aantal geslachten dat op tijdstip t tussen x en x + dx tijdseenheden heeft bestaan dan s es (t−x) dx. De kans dat op tijdstip t een willekeurig gekozen geslacht tussen x en x + dx tijdseenheden heeft bestaan is s es,x dx. Als een op tijdstip t willekeurig gekozen geslacht tussen x en x + dx tijdseenheden heeft bestaan, is de kans dat dit geslacht n soorten bevat volgens formule (15.3) gelijk aan e−r x (1 − e−r x )n−1 voor alle n 1. De kans qn dat 2 Yule beschouwde ook het geval waarin t niet erg groot kan worden verondersteld in vergelijking met de verdubbelingstijd van est. De berekeningen zijn iets ingewikkelder maar de eindresultaten zijn niet erg verschillend. 88 een op tijdstip t willekeurig geslacht n soorten bevat, is dus
  
  Z +∞ qn = s e−s x e−r x (1 − e−r x )n−1 dx. 0
  
  Stel u = r/s. Een eenvoudige rekensom leert dat q1 = 1/(1 + u) en dat qn = 1 u 2u (n − 1)u ··· 1 + u 1 + 2u 1 + 3u 1 + nu (15.5) voor alle n 2. Er geldt immer dat (1 − e−r x )n−1 = (1 − e−r x )n−2 (1 − e−r x ). Dus qn = qn−1 − s Z +∞ e−(r+s) x (1 − e−r x )n−2 e−r x dx. 0 Als we partieel integreren, krijgen we qn = q
  
  n
  
  −1
  
  −
  
  
  
  r+s qn
  
  (
  
  n − 1) r (n − 1) r/s qn−1. 1 + n r/s
  
  en qn = Uit formule (15.5) blijkt dat de rij van kansen (qn )n1 een dalende rij is. Het maximum wordt dus bereikt voor n = 1: de meeste geslachten bevatten slechts één soort. Dit is precies wat de gegevens hadden aangetoond. Bovendien is de afname van qn naar 0 als n naar oneindig gaat, relatief langzaam omdat qn /qn−1 → 1. Dit kan verklaren waarom sommige geslachten een groot aantal soorten bevatten. Om precies te zijn toonde Yule aan dat log qn lineair afneemt met log(n). Z +∞ Introduceer Euler’s Gamma functie Γ(z) = t z−1 e−t dt. Dan Γ(n + 0 1) = n! = n × (n − 1) × · · · × 2 × 1 als n een geheel getal is en Γ(z + 1) = z Γ(z). Dus (15.5) heeft de vorm qn = (n − 1)! u(1 + 1u )(2 + 1u ) · · · (n + u1 ) = Γ(n) Γ(1 + 1u ) u Γ(n + 1 + u1 ). Maar Stirling’s benadering geeft log Γ(n) ≈ n log n−n− 21 log n+constante. Op dezelfde manier, log Γ(n + 1 + 1/u) ≈ n log n − n + ( 1u + 12 ) log n + constante. Ten slotte log qn ≈ −(1 + 1u ) log n + constante. Hoofdstuk 15 89 Beschouw bijvoorbeeld het geval van de hagedissen. De parameter u kan geschat worden uit het aandeel q1 = 1/(1 + u) van geslachten die slechts één soort bevatten. Volgens tabel 15.1 hebben we q1 = 105/259 dus u ≈ 1,467. We kunnen dan de theoretische kans qn berekenen; vervolgens kunnen we het verwachte aantal Qn van geslachten die n soorten bevatten, berekenen door qn te vermenigvuldigen met het totale aantal soorten, dat 259 is (Tabel 15.2). Yule merkte op dat de overeenstemming tussen de waarnemingen en de berekeningen relatief goed3 is gezien de eenvoud van het model, dat bijvoorbeeld geen rekening houdt met de cataclysmen die soorten in de loop van miljoenen jaren evolutie hebben doorgemaakt.
  
  Tabel 15.2: Vergelijking tussen gegevens en theorie in het geval van hagedissen (1.580 soorten ingedeeld in 259 geslachten). A
  
  antal soorten per geslacht Waargenomen aantal geslachten Berekend aantal geslachten 1 2 3 4 5 6 7 8 9 10 11-20 21-30 31-40 41-50 51-60 61-70 71-80 81-90 91-100 101- 105 44 23 14 12 7 6 4 5 5 17 9 3 2 0 1 0 0 0 2 105,0 39,2 21,3 13,6 9,6 7,2 5,6 4,5 3,7 3,1 16,6 6,9 3,9 2,6 1,9 1,4 1,1 0,9 0,7 10,1 totaal 259 259 3
  
  Voor het aantal geslachten met meer dan 100 soorten kreeg Yule een betere fit dan in tabel 15.2 door ervan uit te gaan dat t niet groot was in vergelijking met de verdubbelingstijd van est. 90 Na 1931 trok Yule zich geleidelijk terug van de Universiteit van Cambridge. Hij raakte geïnteresseerd in de statistische verdeling van de lengte van zinnen om schrijvers van boeken te identificeren. Hij paste dit in het bijzonder toe op het boek van John Graunt (zie hoofdstuk 2) dat mogelijk geïnspireerd is door William Petty. In 1944 publiceerde hij een boek over De statistische studie van literaire woordenschat. Hij overleed in 1951. Tegenwoordig wordt het model van Yule nog steeds gebruikt voor de analyse van „fylogenetische bomen” (de genealogische bomen van soorten). Deze bomen, die lijken op die in figuur 15.3, zijn beter bekend dankzij de nieuwe gegevens uit de moleculaire biologie. Maar de toepassingen van het stochastisch proces gedefinieerd door vergelijkingen (15.1)-(15.2) zijn niet beperkt tot de evolutietheorie. Dit proces is een bouwsteen van vele modellen in de populatiedynamica, van het microscopisch niveau (om bijvoorbeeld kolonies van bacteriën te modelleren) tot het macroscopisch niveau (om het begin van een epidemie te modelleren). Het wordt wel „zuiver geboorteproces” of „Yule proces” genoemd. Een eenvoudige variant omvat een kans m dt op overlijden gedurende een klein tijdsinterval dt: de verwachte populatiegrootte op tijdstip t voor dit „geboorte- en sterfteproces” is dan e(r−m)t. De kansverdeling (15.5) wordt ook wel de Yule-verdeling genoemd. Verdelingen met staarten die voldoen aan machtswetten hebben veel aandacht getrokken in verschillende gebieden van de wetenschap. De studie naar epidemieën in stochastische netwerken waar de graden verveeld zijn volgens een machtswet, is slechts één voorbeeld. Verdere
  
  lectuur 1. Aldous, D.J.: Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today. Stat. Sci. 16, 23–34 (2001) projecteuclid.org 2. Edwards, A.W.F.: George Udny Yule. In: Heyde, C.C., Seneta, E. (eds.) Statisticians of the Centuries, 292–294. Springer (2001) 3. McKendrick, A.G.: Studies on the theory of continuous probabilities with special reference to its bearing on natural phenomena of a progressive nature. Proc. Lond. Math. Soc. 13, 401–416 (1914) 4. Simon, H.A.: On a class of skew distribution functions. Biometrika 42, 425– 440 (1955) 5. Willis, J.C.: Age and Area. Cambridge (1922) archive.org 6. Yates, F.: George Udny Yule. Obit. Not. Fellows R. Soc. 8, 308–323 (1952) 7. Yule, G.U.: A mathematical theory of evolution, based on the conclusions of Dr. J. C. Willis, Phil. Trans. Roy. Soc. Lond. B 213, 21–87 (1925) gallica.bnf.fr
  
  Hoofdstuk 16 McKendrick en Kermack over epidemische modellering (1926–1927)
  
  In 1926 bestudeerde McKendrick een stochastisch epidemiemodel en vond een methode om de kans te berekenen dat een epidemie een bepaalde uiteindelijke omvang bereikt. Hij ontdekte ook de partiële differentiaalvergelijking voor leeftijdsgestructureerde populaties in een continue-tijdskader. In 1927 bestudeerden Kermack en McKendrick een deterministisch epidemiemodel en verkregen zij een vergelijking voor de uiteindelijke omvang van de epidemie, waaruit blijkt dat er een bepaalde drempelwaarde is voor de bevolkingsdichtheid. Grote epidemieën kunnen wel boven, maar niet onder deze drempel optreden. Deze ideeën worden nog steeds veel gebruikt in de hedendaagse epidemiologie. Anderson Gray McKendrick werd in 1876 in Edinburgh geboren als laatste van vijf kinderen. Hij studeerde geneeskunde aan de universiteit van Glasgow waar zijn vader professor in de fysiologie was. In 1900 trad hij in dienst bij de Indiase Medische Dienst. Voordat hij naar India ging, vergezelde hij Ronald Ross op een missie ter bestrijding van malaria in Sierra Leone. Daarna diende hij 18 maanden in het leger in Soedan. Bij zijn aankomst in India werd hij aangesteld als arts in een gevangenis in Bengalen waar hij dysenterie onder controle probeerde te krijgen. In 1905 trad hij in dienst bij het nieuwe Centraal Instituut voor Medisch Onderzoek in Kasauli (in het noorden van India). Hij werkte aan hondsdolheid maar studeerde ook wiskunde. In 1920 keerde hij na besmetting met een tropische ziekte terug naar Edinburgh en werd hij hoofd van het laboratorium van het Royal College of Physicians. In 1926 publiceerde McKendrick een artikel over de Toepassingen van de wiskunde op medische problemen, dat verscheidene nieuwe ideeën bevatte. Hij introduceerde met name een continue-tijd wiskundig model voor epidemieën dat rekening hield met het stochastische aspect van infectie en herstel. Beschouw een populatie van grootte N met aanvankelijk slechts één besmette persoon. Mensen kunnen achtereenvolgens drie toestanden doorlopen: de vatbare toestand S, de geïnfecteerde toestand I en de herstelde toestand R 92
  
  Figuur 16.1: McKendrick (1876–1943) en Kermack (1898–1970) (figuur 16.2)1. S I R Figuur 16.2: Mogelijke staten: vatbaar (S), geïnfecteerd (I), hersteld (R). Zij pi,r (t) de kans dat de bevolking op tijdstip t precies i personen in toestand I en r personen in toestand R bevat, waarbij i en r gehele getallen zijn zodanig dat 1 i + r N. In dit geval wordt gezegd dat de bevolking zich in toestand (i, r) bevindt. Het aantal vatbare personen is dan s = N − i − r. In navolging van het werk van Ross over malaria (zie hoofdstuk 12), veronderstelde McKendrick dat gedurende een klein tijdsinterval dt de kans op een nieuwe infectie gelijk is aan a s i dt (d.w.z. evenredig met zowel het aantal vatbare personen als het aantal geïnfecteerde personen). De kans op één nieuwe genezing is gelijk aan b i dt. Zowel a als b zijn positieve parameters. Om pi,r (t + dt) te berekenen, moeten verschillende gevallen worden onderscheiden: • de populatie bevindt zich in toestand (i − 1, r) op tijdstip t en een nieuwe infectie brengt de populatie in toestand (i, r) tussen t en t + dt; 1 Het model van Daniel Bernoulli (zie hoofdstuk 4) omvatte de toestanden S en R maar niet I, omdat de duur van de besmetting veel korter is dan de gemiddelde levensverwachting. de waarschijnlijkheid van deze gebeurtenis is a s (i − 1) dt met s = N − (i − 1) − r; • de populatie bevindt zich in toestand (i, r) op tijdstip t en een nieuwe infectie brengt de populatie in toestand (i + 1, r) tussen t en t + dt; de waarschijnlijkheid van deze gebeurtenis is a s i dt met s = N − i − r; • de populatie bevindt zich in toestand (i + 1, r − 1) op tijdstip t en een nieuw herstel brengt de populatie in toestand (i, r) tussen t en t + dt; de waarschijnlijkheid van deze gebeurtenis is b (i + 1) dt; • de populatie bevindt zich in toestand (i, r) op tijdstip t en één nieuw herstel brengt de populatie in toestand (i − 1, r + 1) tussen t en t + dt; de waarschijnlijkheid van deze gebeurtenis is b i dt. Vandaar dat McKendrick de vergelijkingen d pi,r = a (N − i − r + 1) (i − 1) pi−1,r − a (N − i − r) i pi,r dt + b (i + 1) pi+1,r−1 − b i pi,r (16.1) verkreeg voor 1 i + r N. De eerste term aan het rechterlid ontbreekt als i = 0, terwijl de derde term ontbreekt als r = 0. De beginvoorwaarden zijn pi,r (0) = 0 voor alle (i, r) behalve p1,0 (0) = 1. Met dit model is McKendrick erin geslaagd de kans te berekenen dat de epidemie eindigt met n besmette personen, wat de limiet is van p0,n (t) wanneer t → +∞. Het is inderdaad niet nodig om systeem (16.1) op te lossen. Het volstaat op te merken dat zolang er i geïnfecteerden en r herstelden zijn, de kans op een nieuwe infectie gedurende een klein tijdsinterval dt a (N − i − r) i dt is en de kans op een nieuwe genezing b i dt is. De overgangskansen (zoals ze gewoonlijk genoemd worden in de theorie van Markovketens) van toestand (i, r) naar toestand (i + 1, r) of toestand (i − 1, r + 1) zijn dus respectievelijk P(i,r)→(i+1,r) = a (N − i − r), a (N − i − r) + b P(i,r)→(i−1,r+1) = b, a (N − i − r) + b voor alle i 1 (figuur 16.3). Zij qi,r de kans dat de populatie tijdens de epidemie toestand (i, r) doorloopt. Aangezien i = 1 en r = 0 wanneer t = 0, hebben we q1,0 = 1. De andere toestanden worden bereikt na een besmetting of na een herstel: qi,r =
  
  qi
  
  −1,
  
  r P(i−1
  
  ,
  
  r
  
  )→
  
  (i,r) + qi+1,r−1 P(i+1,r−1)→(i,r).
  
  Figuur 16.3: Diagram met de mogelijke toestanden van een populatie met N = 5 (i op de horizontale as, r op de verticale as) en de mogelijke overgangen ten gevolge van infectie (horizontale pijlen) of van herstel (andere pijlen). i De eerste term van het rechterlid ontbreekt als i = 0 of i = 1. De tweede term ontbreekt als r = 0. Uit deze formule kunnen we eerst (qi,0 )2iN berekenen, dan (qi,1 )0iN−1, dan (qi,2 )0iN−2 enz. De kans dat de epidemie uiteindelijk n mensen zal besmetten is q0,n. In 1926 waren dergelijke berekeningen zeer omslachtig. Daarom beperkte McKendrick zich tot voorbeelden van zeer kleine populaties, bijvoorbeeld een gezin. Met N = 5 personen en b/a = 2 verkreeg hij Tabel 16.1. De grootste waarschijnlijkheden komen overeen met het geval waarin slechts één persoon in het gezin besmet is en met het geval waarin het hele gezin besmet is.
  
  Tabel 16.1: Waarschijnlijkheid dat een epidemie in een gezin van vijf personen n besmet wanneer b/a = 2. n 1 2 3 4 5 q0,n 0,33 0,11 0,09 0,13 0,34
  
  Hetzelfde artikel uit 1926 bevat ook een nieuwe formulering van demografische problemen wanneer tijd als een continue variabele wordt beschouwd. Voor dx oneindig klein, zij P(x,t)dx de bevolking met een leeftijd tussen x en x + dx op tijdstip t. Zij m(x) het sterftecijfer op leeftijd x. Dan is P(x + h,t + h) ≈ P(x,
  
  t
  
  ) −
  
  m(x) P(x,t) h Hoofd
  
  stuk
  
  16 95
  
  voor h oneindig klein. We introduceren de partiële afgeleiden van de functie P(x,t):
  
  ∂P P(x + h,t) − P(x,t) (x,t) = lim, h→0 ∂x h ∂P P(x,t + h) − P(x,t) (x,t) = lim. h→0 ∂t h Met behulp van de ben
  
  adering P(x + h,t + h) ≈ P(x,t) + h ∂P ∂P (x,t) + h (x,t), ∂x ∂t verkreeg McKendrick de volgende partiële differentiaalvergelijking: ∂P ∂P (x,t) + (x,t) + m(x) P(x,t) = 0. ∂t ∂
  
  x Een dergelijke vergelijking komt ook voor bij bevolkingsproblemen die gestructureerd worden door een continue variabele, zoals leeftijd in de demografie (zie hoofdstuk 25) of tijd sinds besmetting in de epidemiologie. In 1921 werd William Ogilvy Kermack benoemd tot hoofd van de chemische afdeling van het Royal College of Physicians Laboratory in Edinburgh. Kermack werd in 1898 geboren in een klein stadje in Schotland. Hij studeerde aan de Universiteit van Aberdeen en begon onderzoek te doen op het gebied van de organische scheikunde in een industrieel laboratorium in Oxford. Ondanks het feit dat hij volledig blind werd na een explosie in zijn laboratorium in Edinburgh in 1924, zette hij zijn chemisch werk voort met de hulp van collega’s en studenten. Kermack begon ook samen te werken met McKendrick aan de wiskundige modellering van epidemieën. Vanaf 1927 publiceerden zij samen een reeks Bijdragen tot de mathematische theorie van epidemieën waarin zij deterministische epidemiemodellen bestudeerden. Stel N is de bevolkingsgrootte met N groot genoeg. Veronderstel zoals in het artikel van 1926 dat mensen ofwel vatbaar, besmet of hersteld kunnen zijn. Indien de ziekte dodelijk is, is de derde toestand in feite de dood. Stel S(t), I(t) en R(t) zijn het aantal mensen in elk van de drie toestanden. Het model is (in vereenvoudigde vorm) een stelsel van drie differentiaalvergelijkingen: dS = −a S I, dt dI = a S I − b I, dt dR = b I. Het aantal nieuwe besmettingen per tijdseenheid is dus, zoals in het stochastische model van 1926, evenredig met zowel het aantal vatbare personen als het aantal besmette personen. Bij het begin van de epidemie, op tijdstip t = 0, is een bepaald aantal mensen besmet: S(0) = N − I0, I(0) = I0 en R(0) = 0, uitgaande van 0 < I0 < N. Hoewel het systeem (16.2)-(16.4) geen gesloten oplossing heeft, kunnen de volgende conclusies getrokken worden: • de totale bevolking S(t) + I(t) + R(t) blijft constant en is gelijk aan N; • S(t), I(t) en R(t) blijven nonnegatief (zoals het hoort aangezien het om populaties gaat); • wanneer t → +∞, S(t) afneemt tot een limiet S∞ > 0, I(t) naar 0 neigt en R(t) toeneemt tot een limiet R∞ < N; • en bovendien de formule − log S∞ a = (N − S∞ ), S(0) b (16.5) geeft impliciet S∞ en dus ook de uiteindelijke epidemieomvang R∞ = N − S∞. Inderdaad, we zien eerst dat dtd (S + I + R) = 0. Dus S(t) + I(t) + R(t) = S(0) + I(0) + R(0) = N. Vergelijkingen (16.2) en (16.3) kunnen worden herschreven als i Rt dh S(t) ea 0 I(τ) dτ = 0, dt i Rt dh I(t) ebt−a 0 S(τ) d
  
  τ
  
  = 0.
  
  dt
  
  Hieruit volgt dat
  
  Rt S(t)
  
  
  
  = S(0) e−a 0 I(τ) dτ > 0, Rt I(t) = I(0) ea 0 S(τ) dτ−bt > 0. Uit vergelijkingen (16.2) en (16.4) blijkt dan dat de functie S(t) dalend is en dat de functie R(t) stijgend is (in het bijzonder, R(t) 0). Daar S(t) 0 en R(t) N, hebben de functies S(t) en R(t) wel degelijk limieten als t → +∞. Daar I(t) = N − S(t) − R(t), heeft I(t) ook een limiet als t → +∞, die alleen nul kan zijn zoals blijkt uit de integratie van (16.4). Uit vergelijking (16.2) blijkt ook dat − d [log S] = a I. dt Als we integreren tussen t = 0 en t = +∞, vinden we log S(0) − log S∞ = a Z +∞ I(t) dt. 0 Vergelijking (16.3) kan worden herschreven als dI dS = − − b I. dt dt Als we integreren tussen t = 0 en t = +∞, krijgen we −I(0) = S(0) − S∞ − b Z +∞ I(t) dt. 0 Door de twee resultaten te combineren bekomen we formule (16.5), waaruit blijkt dat S∞ > 0. Wanneer het aanvankelijke aantal geïnfecteerden I0 klein is in verhouding tot de bevolkingsomvang N, wat vaak het geval is bij het begin van een epidemie in een stad, kan formule (16.5) met S∞ = N − R∞ worden herschreven als R∞ R∞ − log 1 − ≈ R0, (16.6) N N waarbij per definitie aN. b Vergelijking (16.6) heeft alleen een positieve oplossing als R0 > 1. Kermack en McKendrick komen dus tot de volgende conclusie: de epidemie infecteert slechts een niet-verwaarloosbaar deel van de bevolking als R0 > 1. Er is een drempel voor de bevolkingsdichtheid N ∗ = b/a waaronder geen epidemieën kunnen voorkomen. Wanneer de populatiegrootte N net boven deze drempel ligt (N = N ∗ + ε), treedt een epidemie van kleine amplitude op. Uit (16.6) volgt dat R∞ ≈ 2 ε. Dus S∞ ≈ N ∗ − ε: de epidemie brengt de vatbare populatie even ver onder de drempel N ∗ als ze er aanvankelijk boven lag. R0 = 98 2 Als we de benadering − log(1 − x) ≈ x + x2 gebruiken, wordt vergelijking (16.6) R∞ R∞ 1 R∞ 2 + ≈ R0. N 2
  
  N N Dus R∞ ≈ 2(R0 − 1)N = 2 ε (N ∗ + ε) ≈ 2ε. N∗
  
  Zoals in het malariamodel van Ross (hoofdstuk 12) heeft de voorwaarde R0 > 1 een eenvoudige interpretatie. Aangezien a N het aantal personen is dat één besmette persoon per tijdseenheid besmet bij het begin van de epidemie en aangezien 1/b de gemiddelde besmettelijke periode is, is R0 = a N/b het gemiddeld aantal secundaire gevallen ten gevolge van één besmette persoon bij het begin van de epidemie. Voor dodelijke ziekten is R(t) het cumulatieve aantal sterfgevallen sinds het begin van de epidemie en dR/dt het aantal sterfgevallen per tijdseenheid. Kermack en McKendrick hebben vastgesteld dat de grafiek van de functie dR/dt in hun wiskundig model de klokvorm heeft die men van een epidemiecurve verwacht (figuur 16.4). 1 000 doden per week 800 600
  
  Figuur 16.4: De curve dR/dt als functie van de tijd en de gegevens over het aantal sterfgevallen per week tijdens een pestepidemie in Bombay in 1905-1906. 400 200 weken 0 0 5 10 15 20 25 30
  
  Om dR/dt af te leiden deelden zij (16.2) door (16.4) om dS/dR = −a S/b te verkrijgen. Dus S(t) = S(0) exp
  
  (
  
  − R(t)/b). Door dit in vergelijking (16.4) te vervangen en met behulp van S(t) + I(t) + R(t) = N, krijgen ze de vergelijking h a i dR = b N − R − S(0) exp − R, dt b (16.7) die nog steeds niet expliciet kan worden opgelost. Als R(t) echter gedurende de hele epidemie klein blijft, geeft de benadering exp(−u) ≈ 1 − u + u2 /2 dR a a2 ≈ b N − R − S(0) + S(0) R − S(0) 2 R2. (16.8) dt b 2b Dit is een zogenaamde Riccati vergelijking met twee constante oplossingen, een positieve R+ en een negatieve R−, gegeven door de wortels van de tweede-orde polynoom in R aan het rechterlid van (16.8). Zij R̃(t) de exacte oplossing van (16.8) en stel Q(t) = R̃(t) − R+. Dan voldoet Q(t) aan een Bernoulli differentiaalvergelijking gelijk aan die van Daniel Bernoulli en Verhulst (zie (4.5) en (6.1)). Men kan dus formule (6.2) direct aanpassen om Q(t) te krijgen. Een gemakkelijke maar vervelende berekening leert dat dQ/dt van de vorm α cosh2 (βt − γ) is, waarbij α, β en γ constanten zijn die op een ingewikkelde manier afhangen van de parameters van het model. Aangezien dR/dt ≈ R̃/dt = dQ/dt, konden Kermack en McKendrick (α, β, γ) kiezen om hun gegevens te passen. Moderne computers en software kunnen natuurlijk gemakkelijk de differentiaalvergelijking (16.7) numeriek oplossen zonder deze benaderingen te hoeven te doen. De aldus verkregen curve voor dR/dt paste goed bij de gegevens over het aantal sterfgevallen per week tijdens de pestepidemie in Bombay tussen december 1905 en juli 1906 (figuur 16.4) Kermack en McKendrick beschouwden ook het meer algemene model waarbij de besmettelijkheid a(x) afhangt van de tijd x sinds de besmetting en waarbij het genezingspercentage b(x) eveneens van x afhangt. De vergelijking die de uiteindelijke omvang van de epidemie geeft (wanneer het aanvankelijke aantal geïnfecteerde gevallen klein is) is nog steeds (16.6), maar 100 met R0 = N Z +∞ Rx a(x) e− 0 b(y) dy dx. (16.9) 0 De parameter R0 heeft dezelfde interpretatie als in het vorige geval: het is het gemiddeld aantal secundaire gevallen veroorzaakt door één besmette persoon aan het begin van de epidemie. Merk de overeenkomst op tussen (16.9) en de formule van Lotka (10.2) voor R0 in de demografie: leeftijd wordt vervangen Rx door tijd sinds besmetting, overleving door de kans e− 0 b(y) dy om nog besmet te zijn, vruchtbaarheid door de contactfrequentie N a(x). Kermack en McKendrick ontwikkelden in de jaren 1930 verschillende andere wiskundige modellen van epidemieën. Deze vormen nog steeds de bouwstenen voor de meeste van de complexere modellen die tegenwoordig in de epidemiologie worden gebruikt. De parameter R0 speelt nog steeds een centrale rol in de analyse van het model. McKendrick ging in 1941 met pensioen en overleed in 1943. Tussen 1930 en 1933 was Kermack co-auteur van enkele artikelen over mathematische fysica met William McCrea en Edmund Whittaker, beiden van de afdeling wiskunde van de Universiteit van Edinburgh. In de jaren 1930 en 1940 probeerde Kermack’s team van scheikundigen nieuwe moleculen met antimalaria werking te synthetiseren, maar met beperkt succes. In 1938 schreef Kermack samen met Philip Eggleton een populair boek over elementaire biochemie, Het materiaal waar we van gemaakt zijn. In 1944 werd hij tot fellow van de Royal Society gekozen en in 1949 aanvaardde hij de leerstoel biochemie aan de Universiteit van Aberdeen. Later werd hij decaan van de Faculteit der Wetenschappen. Hij ging in 1968 met pensioen en overleed in 1970. Verdere lectuur
  
  1. Advisory Committee appointed by the Secretary of State for India, the Royal Society and the Lister Institute: Reports on plague investigations in India, XXII. J. Hyg. 7, 724–798 (1907) ncbi.nlm.nih.gov 2. Davidson, J.N., Yates, F., McCrea, W.H.: William Ogilvy Kermack 1898– 1970. Biog. Mem. Fellows R. Soc. 17, 399–429 (1971) 3. Gani, J.: A.G. McKendrick. In: Heyde, C.C., Seneta, E. (eds.) Statisticians of the Centuries, 323–327. Springer, New York (2001) 4. Harvey, W.F.: A.G. McKendrick 1876–1943. Edinb. Med. J. 50, 500–506 (1943) 5. McKendrick, A.G.: Applications of mathematics to medical problems. Proc. Edinb. Math. Soc. 13, 98–130 (1926) 6. Kermack, W.O., McKendrick, A.G.: A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. A 115, 700–721 (1927) gallica.bnf.fr
  
  Hoofdstuk 17 Haldane en mutaties (1927)
  
  In een ander deel van zijn artikel uit 1922 beschouwde Fisher het probleem van een mutant gen dat kan worden overgedragen op een willekeurig aantal nakomelingen met een gegeven waarschijnlijkheidsverdeling. Het probleem was formeel hetzelfde als dat van het uitsterven van familienamen, maar dan in een genetische context. Fisher toonde aan dat als de kansverdeling een Poisson-verdeling was en als het mutante gen geen selectief voordeel had, het mutante gen zeer langzaam uit de populatie kon verdwijnen. In 1927 zette de Britse bioloog Haldane de studie van dit model verder en toonde aan dat de waarschijnlijkheid dat een mutant gen zich zou handhaven tweemaal zijn selectief voordeel was. Hij gaf ook een meer rigoureuze behandeling van het uitstervingsprobleem. John Burdon Sanderson Haldane werd in 1892 geboren in Oxford, waar zijn vader professor in de fysiologie was aan de universiteit. Haldane studeerde aan het Eton College en na 1911 aan het New College van de Universiteit van Oxford. Nadat hij zich in zijn eerste jaar op wiskunde had toegelegd, richtte hij zich op de menswetenschappen. Zijn studies werden onderbroken door de Eerste Wereldoorlog, waarin hij diende in Frankrijk en Irak. Nadat hij gewond was geraakt, werd hij als militair instructeur naar India gestuurd. In 1915 publiceerde hij een eerste artikel over genetische experimenten op muizen waarmee hij voor de oorlog was begonnen. In 1919 werd hij fellow van New College, waar hij fysiologie doceerde en net als zijn vader de ademhaling bestudeerde. In 1923 trad hij toe tot het biochemisch laboratorium van F. G. Hopkins1 aan de universiteit van Cambridge, waar hij zich toelegde op de kinetica van enzymen. Hij publiceerde ook een science fiction roman, Daedalus of Wetenschap en de Toekomst (1923), en een essay getiteld Calinicus, Een verdediging van chemische oorlogsvoering (1925). Tussen 1924 en 1934 schreef hij een serie van tien artikelen onder de titel Een wiskundige theorie van natuurlijke en kunstmatige selectie. 1 Frederick Gowland Hopkins, die in 1929 de Nobelprijs voor fysiologie of geneeskunde kreeg voor zijn werk over vitamines. 102
  
  Figuur 17.1: Haldane (1892–1964)
  
  In het vijfde artikel van de reeks, gepubliceerd in 1927, heroverwoog Haldane een ander genetisch model dat Fisher in 1922 had bestudeerd, een model gericht op mutaties. Fisher had de waarschijnlijkheid bestudeerd dat een mutant gen een populatie invadeert of verdwijnt. Dit probleem is formeel hetzelfde als dat van Bienaymé, Galton en Watson betreffende het uitsterven van familienamen. Maar Fisher verwees niet naar deze werken, hoewel hij misschien het artikel van Galton en Watson heeft gelezen dat in de appendix van Galton’s boek Natuurlijke erfenis uit 1889 staat. Zoals in hoofdstuk 9, noem pk de waarschijnlijkheid dat een gen wordt overgedragen op k nakomelingen in de eerste generatie (k 0). Fisher beschouwde ook de genererende functie f (x) = p0 + p1 x + p2 x2 + · · · + pk xk + · · ·, behalve dat hij geen bovengrens voor k vaststelde: de som kan een oneindig aantal termen bevatten. Hij besefte dat, uitgaande van één individu met het mutante gen in generatie 0, de kans dat dit gen zich in k individuen bevindt, de coëfficiënt is van xk in f1 (x) = f (x) voor generatie 1, in f2 (x) = f ( f (x)) voor generatie 2, in f3 (x) = f ( f ( f (x))) voor generatie 3 enz. Op deze manier wordt duidelijk dat de volgende vergelijking geldt. fn (x) = f ( fn−1 (x)) (17.1) Deze vergelijking is veel praktischer dan de door Watson afgeleide vergelijking fn (x) = fn−1 ( f (x)). In het bijzonder volgt uit (17.1) dat de uitsterfkans binnen n generaties xn = fn (0) voldoet aan de iteratieformule xn = f (xn−1 ), zoals Bienaymé reeds had opgemerkt. Als voorbeeld beschouwt Fisher het geval van een plant met een gemuteerd gen dat N zaden kan voortbrengen, waarbij elk zaadje een waarschijnlijkheid q heeft om te overleven en een nieuwe plant voort te brengen. De Hoofdstuk 17 103 kans pk op k nakomelingen met het gemuteerde gen is binomiaal: N k q (1 − q)N−k pk = k voor alle 0 k N en pk = 0 voor k > N. De genererende functie is dan f (x) = (1 − q + q x)N. Laat R0 = N q het gemiddelde aantal zaden zijn dat overleeft om een nieuwe plant voort te brengen. Als N groot is en q klein, dan N +∞ R0 (R0 x)k f (x) = 1 + (x − 1). ≈ eR0 (x−1) = e−R0 ∑ N k! k=0
  
  De kansverdeling (pk ) neigt naar e−R0 (R0 )k /k!, wat een Poisson-verdeling wordt genoemd. Fisher berekende dan de uitstervingskans binnen n generaties, gebruik makend van x0 = 0, xn ≈ eR0 (xn−1 −1) en de numerieke waarden N = 80 en q = 1/80. In dit geval is R0 = Nq = 1. Een omslachtige berekening leert dat x100 ≈ 0,98: een mutant gen zonder selectief voordeel (R0 = 1) verdwijnt heel langzaam. Er is nog 2 % kans dat het gen na 100 generaties nog in de populatie aanwezig is. In 1922 heeft Fisher de studie van dit model niet voortgezet. Voortbordurend op het werk van Fisher, merkte Haldane in zijn artikel van 1927 voor het eerst op dat voor elke kansverdeling (pk ) zodanig dat p0 > 0, de vergelijking x = f (x) precies twee wortels heeft in het interval (0, 1] wanneer het gemiddelde aantal nakomelingen dat drager is van het mutante gen R0 strikt groter is dan 1, d.w.z. wanneer het mutante gen een selectief voordeel heeft. Bovendien is de uitstervingskans x∞, die de limiet is van xn als n → +∞, de kleinste van de twee wortels van x = f (x): het gen heeft een positieve kans om zich in de populatie te vestigen. In tegenstelling tot Bienaymé en Cournot leverde Haldane een bewijs voor deze conclusie. Inderdaad, f 0 (x) 0 en f 00 (x) 0 op het interval [0, 1]. Met andere woorden, de functie f (x) is niet-degressief en convex. De veronderstellingen f (0) = p0 > 0 en f 0 (1) = R0 = p1 + 2p2 + 3p3 + · · · > 1 impliceren dat de vergelijking f (x) = x precies twee oplossingen heeft in het interval (0, 1]: x = 1 en x∗ zo dat 0 < x∗ < 1. Haldane verwees vervolgens naar een artikel van Gabriel Koenigs uit 1883, waarin werd aangetoond dat als xn = f (xn−1 ) en xn → x∞, dan x∞ = f (x∞ ) en | f 0 (x∞ )| 1. Als f 0 (1) > 1, is de enige mogelijkheid dat x∞ = x∗. Voor het geval van een Poisson-verdeling met f (x) = eR0 (x−1) en R0 net iets groter dan 1, ligt de uitdovingskans x∞ heel dicht bij 1. De vergelijking 104 f (x∞ ) = x∞ is dan equivalent met R0 (x∞ − 1) = log x∞ ≈ (x∞ − 1) − (x∞ − 1)2. 2 Hieruit volgt dat 1 − x∞ ≈ 2(R0 − 1). Haldane concludeerde dat de kans dat het gemuteerde gen niet zal uitsterven twee keer zo groot is als het selectieve voordeel R0 − 1. Zonder Haldane te citeren nam Fisher in zijn boek uit 1930 als voorbeeld het geval waarin R0 = 1,01, wat een kans van 2 % oplevert dat het mutante gen niet zal uitsterven. Haldane werd in 1932 lid van de Royal Society. Hij verliet Cambridge om professor in de genetica en later biometrie te worden aan het University College in Londen. Hij was toen vooral geïnteresseerd in menselijke genetica: schatting van mutatiepercentages, genetische kaarten van chromosomen enz. Naast zijn wetenschappelijke boeken (Dierlijke biologie in 1927 met Julian Huxley, Enzymen in 1930 en De oorzaken van de evolutie in 1932, De biochemie van de genetica in 1954), publiceerde hij een groot aantal artikelen over wetenschap in de pers (o.a, over de oorsprong van het leven) en enkele essays (De ongelijkheid van de mens in 1932, De filosofie van een bioloog in 1935, De marxistische filosofie en de wetenschappen in 1938, Erfelijkheid en Politiek in 1938 en De vooruitgang van de wetenschap in 1947). Na verscheidene bezoeken aan Spanje tijdens de burgeroorlog, probeerde hij zijn eigen land ervan te overtuigen schuilkelders te bouwen tegen luchtbombardementen. Tijdens de Tweede Wereldoorlog werkte hij aan ademhalingsproblemen in onderzeeërs. Lid van de communistische partij sinds 1942, nam hij in 1950 ontslag vanwege de officiële afwijzing van de Mendeliaanse genetica in de USSR door de invloed van Lysenko. In 1957 vestigde hij zich in India, waar hij zijn onderzoek voortzette, eerst aan het Indiase Statistische Instituut in Calcutta en later in Bhubaneswar. Nadat hij Indisch staatsburger was geworden, overleed hij in 1964. Verdere lectuur
  
  1. Clark, R.: J.B.S., The Life and Work of J.B.S. Haldane. London (1968) 2. Haldane, J.B.S.: A mathematical theory of natural and artificial selection, Part V, Selection and mutation. Proc. Camb. Philos. Soc. 23, 838–844 (1927) 3. Haldane, J.B.S.: The Causes of Evolution. Longmans (1932) archive.org 4. Pirie, N.W.: John Burdon Sanderson Haldane 1892-1964. Biog. Mem. Fellows R. Soc. 12, 218–249 (1966)
  
  Hoofdstuk 18 Erlang en Steffensen over het uitstervingsprobleem (1929–1933)
  
  In 1929 boog de Deense telefooningenieur Erlang zich opnieuw over het probleem van het uitsterven van familienamen. Zijn landgenoot, de statisticus Steffensen, werkte een volledige oplossing van het probleem uit. Hij toonde met name aan dat de verwachting van het aantal nakomelingen in elke generatie exponentieel toeneemt en sloeg daarmee de brug tussen stochastische en deterministische populatiemodellen. Agner Krarup Erlang werd in 1878 geboren in Lønborg, Denemarken. Zijn vader was schoolmeester. Tussen 1896 en 1901 studeerde de jonge Erlang wiskunde, natuurkunde en scheikunde aan de Universiteit van Kopenhagen. Daarna gaf hij enkele jaren les op middelbare scholen, terwijl hij zich bleef interesseren voor wiskunde, vooral voor kansrekening. Hij ontmoette Jensen, hoofdingenieur bij de Telefoonmaatschappij van Kopenhagen en een amateur-wiskundige, die hem in 1908 overhaalde om bij het nieuwe onderzoekslaboratorium van de maatschappij te komen werken. Erlang begon artikelen te publiceren over de toepassingen van de waarschijnlijkheidsrekening op het beheer van telefoongesprekken. In 1917 ontdekte hij een formule voor wachttijden, die al snel door telefoonmaatschappijen over de hele wereld werd gebruikt. Zijn artikelen, die eerst in het Deens werden gepubliceerd, werden vervolgens in verschillende andere talen vertaald.