Informatie

Hoe kunnen we de constanten van het SIR-model berekenen?


Ik ben geen bioloog van beroep, maar ik probeer een eenvoudig SIR-model uit het volgende artikel te implementeren. Het papier gaat ervan uit dat de volgende constanten bekend zijn:

rc = groeisnelheid van vatbare groep S

k = draagkracht van de gevoelige groep bij afwezigheid van infectie

alfa = maximale waarde van het reductiepercentage per hoofd van S als gevolg van I

een = halve verzadigingsconstante

gamma = natuurlijk herstelpercentage van infectie

delta = sterftecijfer van herstelde populatie

Ik vond enkele coronavirusgegevens uit de volgende bron, en ik heb met name de gegevens van Italië uitgezet:

  • Kan ik bovenstaande constanten uit dit soort gegevens berekenen?
  • Is het een goede aanname dat deze "variabelen" constant zijn?

Hier is mijn poging van het model met de constanten die in het papier worden gegeven


Hoe kunnen we de constanten van het SIR-model berekenen? - Biologie

Deel 2: Het differentiaalvergelijkingsmodel

Als eerste stap in het modelleringsproces identificeren we de onafhankelijke en afhankelijke variabelen. De onafhankelijke variabele is tijd t, gemeten in dagen. We beschouwen twee gerelateerde sets van afhankelijke variabelen.

De eerste set afhankelijke variabelen telt mensen in elk van de groepen, elk als functie van de tijd:

De tweede set afhankelijke variabelen vertegenwoordigt de fractie van de totale bevolking in elk van de drie categorieën. Dus indien N is de totale populatie (7.900.000 in ons voorbeeld), we hebben

    Onder de aannames die we hebben gemaakt, beschrijf in woorden hoe u denkt NS) moet variëren met de tijd? Hoe moet r(t) variëren met de tijd? Hoe moet het) variëren met de tijd?

Vervolgens maken we enkele veronderstellingen over de veranderingssnelheden van onze afhankelijke variabelen:

Niemand is toegevoegd voor de vatbare groep, aangezien we geboorten en immigratie negeren. De enige manier waarop een individu bladeren de vatbare groep is door besmet te raken. We nemen aan dat de tijd-snelheid van verandering van NS), de nummer van vatbare personen, hangt af van het aantal reeds vatbare personen, het aantal reeds geïnfecteerde personen en de hoeveelheid contact tussen vatbare personen en geïnfecteerden. Stel in het bijzonder dat elk geïnfecteerd individu een vast aantal heeft B van contacten per dag die voldoende zijn om de ziekte te verspreiden. Niet al deze contacten zijn met gevoelige personen. Als we uitgaan van een homogene vermenging van de populatie, fractie van deze contacten die met vatbare personen zijn, is: NS). Dus gemiddeld genereert elk geïnfecteerd individu bs(t) nieuwe geïnfecteerde personen per dag.

Laten we eens kijken wat deze aannames ons vertellen over afgeleiden van onze afhankelijke variabelen.

    De vatbare vergelijking .

(a) Leg zorgvuldig uit hoe elk onderdeel van de differentiaalvergelijking

(1)

volgt uit de tekst die aan deze stap voorafgaat. Vooral,

    Waarom is de factor van Het) Cadeau?

(b) Leg nu uit hoe deze vergelijking leidt tot de volgende differentiaalvergelijking voor NS).

(3)

volgt uit een van de aannames voorafgaand aan stap 4.

(4)

Welke veronderstelling over het model weerspiegelt dit?

(b) Leg nu zorgvuldig uit hoe elk onderdeel van de vergelijking

(5)

volgt uit wat je tot nu toe hebt gedaan. Vooral,

Ten slotte voltooien we ons model door elke differentiaalvergelijking een beginvoorwaarde te geven. Voor dit specifieke virus -- Hong Kong-griep in New York City eind jaren zestig -- was aan het begin van de epidemie bijna niemand immuun, dus bijna iedereen was vatbaar. We gaan ervan uit dat er een spoor van infectie was in de populatie, laten we zeggen 10 mensen. Onze initiële waarden voor de populatievariabelen zijn dus:

S(0) = 7.900.000
ik(0) = 10
R(0) = 0

In termen van de geschaalde variabelen zijn deze beginvoorwaarden:

s(0) = 1
i(0) = 1,27 x 10 - 6
r(0) = 0

(Opmerking: de som van onze startpopulaties is niet precies N , en de som van onze fracties ook niet precies 1. Het infectiespoor is zo klein dat dit geen verschil maakt.) Ons volledige model is

We weten nog geen waarden voor de parameters b en k, maar we kunnen ze schatten en ze zo nodig aanpassen om te passen bij de extra sterftegegevens. We hebben de gemiddelde besmettelijkheidsduur al geschat op drie dagen, dus dat zou k = 1/3 suggereren. Als we veronderstellen dat elke geïnfecteerde om de twee dagen een mogelijk besmet contact zou maken, dan zou b 1/2 zijn. We benadrukken dat dit slechts een gok is. De volgende grafiek toont de oplossingskrommen voor deze keuzes van b en k.

    In stap 1 en 2 heb je je ideeën vastgelegd over hoe de oplossingsfuncties eruit zouden moeten zien. Hoe verhouden die ideeën zich tot de bovenstaande figuur? Vooral,

    Wat vindt u van het relatief lage besmettingsniveau op het hoogtepunt van de epidemie?

Vanwege de standaard keuze van symbolen voor de afhankelijke variabelen wordt dit model het SIR-model genoemd. In deel 3 zullen we zien hoe oplossingskrommen kunnen worden berekend, zelfs zonder formules voor de oplossingsfuncties. [Opmerking: als u de module Systemen van differentiaalvergelijkingen al hebt doorlopen, kunt u deel 3 van deze module overslaan en direct naar deel 4 gaan.]


SIR-model

Terwijl onze modellen worden gemotiveerd door een probleem in de neurowetenschappen en terwijl we verwijzen naar onze modellen N als 'neuronale netwerken' is er niets inherent 'neuronaal' aan deze structuren. Het zijn gewoon wiskundige objecten. In het algemeen kan een bepaald wiskundig object, zoals een dynamisch systeem, als model dienen voor een willekeurig aantal zeer verschillende natuurlijke systemen.

Stel je een infectieziekte voor die zich verspreidt in een populatie van N individuen. Op elk willekeurig moment kan een persoon in een van de drie categorieën vallen: geïnfecteerd en vatbaar voor het infecteren van anderen (in de set l), gezond maar vatbaar voor de infectie (in de set S), of hersteld en (tijdelijk of permanent) immuun voor de ziekte (in de set R). Infectie kan het gevolg zijn van een contact tussen een infectieus en een vatbaar persoon. Deze aannames leiden tot de klassieke SIR-modellen van ziektedynamiek (zie [16] voor een recensie). Merk op dat besmettelijkheid vergelijkbaar is met het afvuren van een neuron: het kan dezelfde toestand op een later tijdstip bij een ander individu veroorzaken. In deze sectie zullen we onderzoeken of onze netwerken geschikte modellen kunnen zijn voor de dynamiek van bepaalde ziekten.

De overeenkomst tussen een wiskundig model en het onderliggende natuurlijke systeem is nooit perfect Wiskundige modellering is altijd gebaseerd op vereenvoudigende aannames. De zogenoemde kunst van wiskundig modelleren is in wezen een talent voor het maken van vereenvoudigende aannames die leiden tot modellen die eenvoudig genoeg zijn om verkenning door computersimulaties of wiskundige methoden mogelijk te maken en toch voldoende details bevatten om realistische voorspellingen te doen over een natuurlijk systeem. De eerste vraag die een modelbouwer moet beantwoorden, is welke aspecten van het natuurlijke systeem het model zou moeten voorspellen. Hier zullen we alleen geïnteresseerd zijn in de vraag hoe het aandeel geïnfecteerde individuen in de populatie in de loop van de tijd verandert.

SIR-modellen zijn er met name in een verscheidenheid aan smaken, er zijn veel details waarmee u rekening moet houden die van ziekte tot ziekte verschillen. Bij daadwerkelijke modellering worden deze details afgeleid uit de beschikbare gegevens en wordt het model geconstrueerd door geschikte aannames uit de gegevens af te leiden. Hier hebben we niet de ruimte om feitelijke gegevens voor een ziekte te beschouwen, maar laten we eens kijken naar enkele modelaannames die men zou kunnen hebben afgeleid:

Een geïnfecteerd persoon zal besmettelijk zijn gedurende een tijdsperiode T 0 met een gemiddelde waarde E ( T 0 ) en een zeer kleine standaarddeviatie.

Het tijdsverloop tussen infectie en infectieus worden is erg klein ten opzichte van E (T 0 ), zodat het verwaarloosbaar kan zijn.

Nadat de infectie niet langer besmettelijk is, blijft een persoon immuun voor de ziekte gedurende een tijdsperiode T 1 met E ( T 1 ) pE ( T 0 ) waarbij P is een positief geheel getal gegeven door de gegevens en de standaarddeviatie σ van T 1 is klein.

Contact tussen een geïnfecteerde en een vatbare persoon tijdens het tijdsinterval van besmettelijkheid zal waarschijnlijk resulteren in een nieuwe infectie Q.

Dit zijn alle veronderstellingen die we hier zullen gebruiken. Maar natuurlijk verdoezelen deze veronderstellingen al veel details die de ziektedynamiek (of misschien niet) significant kunnen beïnvloeden. Voordat u leest hoe we een model bouwen op basis van deze aannames, wilt u misschien een paar minuten de tijd nemen om het volgende te doen:

Oefening 6.16

Noem enkele details die door de bovenstaande veronderstellingen worden genegeerd, maar die een aanzienlijke invloed kunnen hebben op hoe de ziekte zich onder de bevolking zal verspreiden. □

Aangezien de tijdsduur dat een individu besmettelijk blijft een "zeer kleine standaarddeviatie" heeft en het tijdsverloop tussen het moment van overdracht en het worden van infectie verwaarloosbaar kan zijn, kunnen we proberen te werken met een tijddiscreet model waarbij de tijdseenheid wordt gekozen als T0. Beschouw nu een netwerk N = 〈 D , p → , 1 → 〉 met V D = [ n ] en constante p → = ( ​​p , … , p ) . Laten we elke i ∈ [ n ] interpreteren als een individu van een vaste populatie, interpreteer toestand s i ( t ) = 0 als individu l op tijd besmettelijk zijn t, staat s i ( t ) = p als individu l vatbaar zijn voor infectie en 0 < s i ( t ) < p als individu aangeven l immuun zijn voor infectie. Een boog 〈 i , j 〉 ∈ A D geeft aan dat individu l interageert met individuele J. Het is natuurlijk om hier aan te nemen dat AD symmetrisch is, dat wil zeggen dat beide bogen 〈 i , j 〉 , 〈 j , i 〉 zich in AD bevinden of geen van beide is, maar zoals we zullen zien in Project 8 [1] deze veronderstelling is eigenlijk niet nodig.

Zal de dynamiek van N overeenkomen met de dynamiek van de ziekte? Mogelijk. Merk op dat s i ( t ) = 0 s i ( t + τ ) ≠ 0 betekent voor τ ∈ [ p ] , wat mooi overeenkomt met immuun zijn voor de ziekte gedurende een periode van pT 0 .

Echter, in de bovenstaande interpretatie van N Een individu J dat is op tijd vatbaar het zal besmettelijk worden op tijdstip t + 1 wanneer er een boog is 〈 i , j 〉 ∈ A D met individu l op tijd besmettelijk zijn t. Dit zou betekenen dat de transmissiekans Q is 1 en ofwel i , j heeft altijd interactie gedurende een tijdsinterval van lengte T 0 (wanneer 〈 i , j 〉 , 〈 j , i 〉 ∈ AD ), of i , j heeft nooit interactie (wanneer 〈 i , j 〉 , 〈 j , ik 〉 ∉ AD ). Deze aannames zijn te extreem om realistisch te zijn. In de praktijk zal er altijd een element van stochastiek zijn, zowel in de structuur van contacten als in de feitelijke overdracht van de ziekte die het gevolg is van een bepaald contact. We kunnen deze effecten in ons model inbouwen door het als volgt te wijzigen. Neem aan dat gedurende elke tijdsperiode met lengte T 0 twee bepaalde individuen een interactie aangaan met waarschijnlijkheid q 0 , en deze interacties zijn onafhankelijk. Als bij een bepaalde interactie een vatbaar en een besmettelijk individu betrokken is, zal er waarschijnlijk een overdracht plaatsvinden Q. Zo kunnen we de dynamiek van een willekeurig discreet dynamisch systeem N1 modelleren waarvan de bijwerkregels identiek zijn aan het model N hierboven, behalve dat de digraaf D t niet vast zal staan, maar willekeurig opnieuw zal worden getrokken bij elke tijdstap uit de verdeling D n ( q 0 q ) van Erdős-Rényi willekeurige digraphs.

Merk op dat N 1 de verhouding E ( T 1 ) / E ( T 0 ) precies als een geheel getal behandelt en de "kleine standaarddeviatie" van T 1 die we uit de gegevens krijgen als praktisch nul beschouwt. Dit kan worden toegestaan, maar men moet voorzichtig zijn met dergelijke vereenvoudigingen. Project 8 [1] nodigt de lezer uit om te onderzoeken of deze vereenvoudigingen al dan niet tot valse voorspellingen zullen leiden. Meer bepaald stelt het project een suite N 1 , N 2 , N 3 , N 4 van steeds meer uitgebreide en schijnbaar realistischere modellen van hetzelfde natuurlijke systeem. Met name de constructies die voor de modellen N 2 tot en met N 4 worden gebruikt, bevatten de standaarddeviatie van T 1 die wordt genegeerd door model N 1 , en stellen ons in staat om te gaan met situaties waarin de empirisch waargenomen verhouding E ( T 1 ) / E ( T 0 ) is slechts bij benadering, maar niet precies een geheel getal, zoals vrijwel zeker het geval zal zijn voor echte datasets. De meer uitgebreide modellen zijn echter moeilijker te bestuderen. We zullen onderzoeken welke van deze modellen, indien aanwezig, precies het juiste detailniveau bevat.

Zonder het juiste antwoord voor Project 8 [1] te verklappen, laten we er omwille van het argument van uitgaan dat onze simulaties voor model N 3 de voorspellingen van gevestigde modellen uit de literatuur bevestigen. Zou dit betekenen dat model N 3 "goed genoeg" is om nauwkeurige voorspellingen te doen over het beloop van een werkelijke ziekte? Dit soort vragen kan niet met zekerheid bevestigend worden beantwoord. De natuur kan altijd verrassingen bevatten in de vorm van verborgen kenmerken van het echte systeem die de dynamiek beïnvloeden, maar die een modelbouwer misschien over het hoofd heeft gezien. Maar het zou mogelijk zijn om wiskundige stellingen te bewijzen dat een bepaald eenvoudiger model zoals ons discrete model N 3 ons exact hetzelfde antwoord op een bepaalde vraag zal geven dan een meer gedetailleerde, zoals een ODE-versie van de overeenkomstige SIR model. We zullen deze vraag niet behandelen voor ziektemodellen, maar in de volgende sectie zullen we een dergelijke stelling beschrijven voor modellen van bepaalde neuronale netwerken. Stellingen van dit type zouden buitengewoon waardevol zijn, aangezien discrete modellen vaak gemakkelijker te bestuderen zijn, althans door simulaties. In de ziektedynamiek speelt bijvoorbeeld het onderliggende netwerk van contacten een belangrijke rol, waarbij sommige individuen frequent contact hebben met vele anderen, terwijl anderen meer sociaal geïsoleerd kunnen zijn. Dit fenomeen kan eenvoudig in ons raamwerk worden gemodelleerd door te tekenen NS van een gegeven verdeling met ongelijke kansen voor verschillende potentiële bogen 〈 i , j 〉 , maar het is moeilijker om op te nemen in een differentiaalvergelijkingenraamwerk.

We willen er echter op wijzen dat dergelijke stellingen alleen de geschiktheid van een bepaalde vereenvoudiging voor een specifieke reeks vragen over modeldynamiek kunnen aanpakken. Merk op dat er ten minste één aspect van ziektedynamiek is waarin al onze modellen flagrant fout zijn: neem bijvoorbeeld aan dat E ( T 0 ) één week is, wat overeenkomt met één tijdstap in onze discrete modellen. Nu, als je ze te letterlijk behandelt, lijkt het alsof onze modellen zeggen dat elk individu ofwel op maandag ziek wordt en op zondag herstelt, ofwel de hele week gezond blijft. Dit is onzin. Onze modellen worden nutteloos voor ziektedynamiek op tijdschalen van minder dan een week.

Deze laatste observatie brengt ons terug naar de neurowetenschap. Onze modellen zijn geïnspireerd op het fenomeen dynamische clustering waar tijd lijkt te zijn verdeeld in discrete afleveringen van ongeveer gelijke lengte. Sommige neuronen vuren terwijl andere neuronen tijdens een bepaalde episode in rust blijven, en het lidmaatschap van de clusters van vurende neuronen verandert van episode tot episode. Merk op dat dit fenomeen precies lijkt op de dynamiek van een hypothetische ziekte waarbij groepen mensen op maandag ziek worden en op zondag herstellen, waarbij verschillende groepen gedurende verschillende weken ziek zijn en het groepslidmaatschap in de loop van de tijd verandert. Dit soort dingen gebeurt niet met ziekten, maar zoals vermeld in paragraaf 6.2, is het empirisch waargenomen in sommige opnames van echte neuronen. Het is een raadselachtig fenomeen, en men wil natuurlijk weten wat de oorzaak kan zijn.

Een letterlijke lezing van ons discrete model zou dit fenomeen voorspellen, maar dit betekent niet dat ons discrete model verklaart dynamische clustering, omdat dit patroon al is ingebouwd in de aannames van discrete tijdstappen voor onze modellen. Als men wil begrijpen waarom dynamische clustering zal optreden in sommige, maar niet in andere, neuronale netwerken, moet men modellen overwegen die gebaseerd zijn op gekoppelde differentiaalvergelijkingen die geen ingebouwde aanname hebben van discrete episodes en studiecondities op deze modellen waaronder het fenomeen zal optreden . In de volgende sectie worden enkele van dit soort resultaten beschreven en hoe deze zich verhouden tot de bredere context van neurowetenschappen.


Het SIR-model voor de verspreiding van ziekten - het differentiaalvergelijkingsmodel

Als eerste stap in het modelleringsproces identificeren we de onafhankelijke en afhankelijke variabelen. De onafhankelijke variabele is tijd  t,  gemeten in dagen. We beschouwen twee gerelateerde sets van afhankelijke variabelen.

De eerste set afhankelijke variabelen telt mensen in elk van de groepen, elk als functie van de tijd:

S = S(t) is het aantal gevoelig individuen,
ik = I(t) is het aantal besmet individuen, en
R = R(t) is het aantal hersteld individuen.

De tweede set afhankelijke variabelen vertegenwoordigt de fractie van de totale bevolking in elk van de drie categorieën. Dus, als  N  is de totale populatie (7.900.000 in ons voorbeeld), we hebben

s(t) = S(t)/N, het gevoelige deel van de bevolking,
i(t) = I(t)/N, de besmette fractie van de bevolking, en
r(t) = R(t)/N, het herstelde deel van de bevolking.

Het lijkt misschien natuurlijker om met populatietellingen te werken, maar sommige van onze berekeningen zullen eenvoudiger zijn als we in plaats daarvan de breuken gebruiken. De twee sets afhankelijke variabelen zijn evenredig met elkaar, dus beide sets geven ons dezelfde informatie over de voortgang van de epidemie.

    Op basis van de veronderstellingen die we hebben gemaakt, hoe denk je dat  NS)  moet variëren met de tijd? Hoe moet  r(t)  met de tijd variëren? Hoe moet  het)  met de tijd variëren?

Vervolgens maken we enkele veronderstellingen over de veranderingssnelheden van onze afhankelijke variabelen:

Niemand is toegevoegd voor de vatbare groep, aangezien we geboorten en immigratie negeren. De enige manier waarop een individu bladeren de vatbare groep is door besmet te raken. We nemen aan dat de tijd-snelheid van verandering van  NS),  de nummer van de vatbare personen hangt 1 af van het aantal reeds vatbare personen, het aantal reeds geïnfecteerde personen en de hoeveelheid contact tussen vatbare personen en geïnfecteerden. Stel in het bijzonder dat elke geïnfecteerde persoon een vast aantal   . heeftB160 van de contacten per dag die voldoende zijn om de ziekte te verspreiden. Niet al deze contacten zijn met gevoelige personen. Als we uitgaan van een homogene vermenging van de populatie, fractie van deze contacten die gevoelig zijn, is  NS).  Zo genereert elk geïnfecteerd individu gemiddeld  b s(t)  nieuwe geïnfecteerde personen per dag. [Met een grote vatbare populatie en een relatief kleine geïnfecteerde populatie, kunnen we lastige telsituaties negeren, zoals een enkele vatbare die meer dan één geïnfecteerde op een bepaalde dag tegenkomt.]

Laten we eens kijken wat deze aannames ons vertellen over afgeleiden van onze afhankelijke variabelen.

    De vatbare vergelijking . Leg zorgvuldig uit hoe elk onderdeel van de differentiaalvergelijking

(1)

volgt uit de tekst die aan deze stap voorafgaat. Vooral,

    Waarom is de factor  Het)  aanwezig?

(3)

volgt uit een van de aannames voorafgaand aan stap 4.

(4)

Welke veronderstelling over het model weerspiegelt dit? Leg nu zorgvuldig uit hoe elk onderdeel van de vergelijking

(5)

volgt uit wat je tot nu toe hebt gedaan. Vooral,

Ten slotte voltooien we ons model door elke differentiaalvergelijking een beginvoorwaarde te geven. Voor dit specifieke virus -- Hong Kong-griep in New York City eind jaren zestig -- was aan het begin van de epidemie bijna niemand immuun, dus bijna iedereen was vatbaar. We gaan ervan uit dat er een spoor van infectie was in de populatie, laten we zeggen 10 mensen. 2  Onze beginwaarden voor de populatievariabelen zijn dus

S(0) =ه.900.000
I(0) =㺊
R(0) =ـ

In termen van de geschaalde variabelen zijn deze beginvoorwaarden:

s(0) =ف
i(0) =ف.27 x㺊 - 6
r(0) =ـ

(Opmerking: de som van onze startpopulaties is niet precies  N,  en de som van onze breuken is ook niet precies  1.  Het sporenniveau van infectie is zo klein dat dit geen verschil zal maken.) Ons complete model is

We kennen geen waarden voor de parameters  B  en   'Nog niet, maar we kunnen ze schatten en ze zo nodig aanpassen aan de overtollige sterftegegevens. We hebben de gemiddelde besmettelijkheidsperiode al geschat op drie dagen, dus dat zou suggereren  k =ف/3.  Als we raden dat elke geïnfecteerde om de twee dagen een mogelijk besmet contact zou maken, dan  B  zou   zijn1/2.  We benadrukken dat dit slechts een gok is. De volgende grafiek toont de oplossingskrommen voor deze keuzes van "160"B  en  k

    In stap 1 en 2 heb je je ideeën vastgelegd over hoe de oplossingsfuncties eruit zouden moeten zien. Hoe verhouden die ideeën zich tot de bovenstaande figuur? Vooral,

    Wat vindt u van het relatief lage besmettingsniveau op het hoogtepunt van de epidemie?

In deel 3 zullen we zien hoe oplossingskrommen kunnen worden berekend, zelfs zonder formules voor de oplossingsfuncties.

1 Merk op dat we het bijvoeglijk naamwoord "vatbaar" in een zelfstandig naamwoord hebben veranderd. Het is gebruikelijk in de epidemiologie om te verwijzen naar 'vatbare personen', 'geïnfecteerden' en 'herstelden' in plaats van altijd langere zinnen te gebruiken, zoals 'bevolking van vatbare mensen' of zelfs 'de vatbare groep'.

2 'Terwijl' ik(0)'is normaal gesproken klein ten opzichte van' N, we moeten hebben I(0) >ـ om zich een epidemie te laten ontwikkelen. Vergelijking (5) zegt redelijkerwijs dat als ik =ـ op tijd 0 (of wanneer dan ook), dan dI/dt =ـ ook, en er kan nooit enige toename zijn van de 0 niveau van infectie.

David Smith en Lang Moore, "Het SIR-model voor de verspreiding van ziekten - het differentiaalvergelijkingsmodel", Convergentie (december 2004)


3. Methodologie

De gegevens in [29] voor China, Zuid-Korea, India, Australië, de VS, Italië en de staat Texas (gemeenschappen) zijn georganiseerd in de vorm van tijdreeksen waarbij de rijen tijdopnames zijn (van januari tot juni 2020 ), en de drie kolommen zijn, het totale aantal gevallen dat ik totd heb (eerste kolom), het aantal geïnfecteerde personen ID kaart (tweede kolom) en sterfgevallen D d (derde kolom). Bijgevolg kan het aantal verwijderingen worden geschat op basis van de gegevens door:

Omdat we de numerieke oplossingen van ons voorgestelde SIR-model (1) willen aanpassen aan de geregistreerde gegevens van [29] , beschouwen we voor elke dataset (gemeenschap) de beginvoorwaarden in het interval [0,1] en schalen ze door een schaal factor F om de geregistreerde gegevens door visuele inspectie te passen. In het bijzonder zijn de beginvoorwaarden voor de drie populaties zo ingesteld dat S ( 0 ) = 1 (dwz alle individuen worden aanvankelijk als vatbaar beschouwd), I ( 0 ) = R m ( 0 ) = I maxd / f < 1 , waarbij I maxd het maximale aantal geïnfecteerde personen is ID kaart . Bijgevolg zijn de parameters a, b, f en I m a x d worden op basis van trial-and-error en visuele inspecties handmatig aangepast om zo goed mogelijk bij de vastgelegde gegevens te passen. Een voorlopige analyse met behulp van niet-lineaire fittingen om het model aan te passen aan de gepubliceerde gegevens [29] leverde op zijn best inferieure resultaten op dan die verkregen in dit artikel met behulp van onze trial-and-error-aanpak met visuele inspecties, in de zin dat de modeloplossingen dat deden niet zo dicht bij de gepubliceerde gegevens volgen, wat onze benadering in de krant rechtvaardigt. Een belangrijke reden hiervoor is dat de gepubliceerde gegevens (inclusief die in [29] die we hier gebruiken) gegevens zijn uit verschillende landen die verschillende methoden volgen om ze vast te leggen, waarbij niet alle geïnfecteerde personen of sterfgevallen zijn verantwoord.

In deze context, S, ik en Rm ≥਀ bij elk t ≥਀. Systeem (1) kan numeriek worden opgelost om te bepalen hoe de geschaalde (by F) gevoelig S, besmet l en verwijderd Rm populaties (wat we modeloplossingen noemen) evolueren met de tijd, in goede overeenstemming met de geregistreerde gegevens. In het bijzonder, omdat dit systeem eenvoudig is met goed opgevoede oplossingen, hebben we de eerste-orde Euler-integratiemethode gebruikt om het numeriek op te lossen, en een tijdstap h = 200/5000 = 0,04 die overeenkomt met een uiteindelijke integratietijd tF van 200 dagen sinds januari 2020. Dit komt neer op een verdubbeling van het tijdsinterval in de geregistreerde gegevens in [29] en maakt voorspellingen mogelijk tot 100 dagen na januari 2020.

Wat natuurlijk belangrijk is bij het bestuderen van de verspreiding van een virus, is het aantal doden NS en terugvorderingen R op tijd. Omdat deze cijfers niet rechtstreeks door het SIR-model worden geleverd (1), hebben we ze geschat door eerst de gegevens voor sterfgevallen in kaart te brengen D d vs de verwijderingen R m d , waarbij R m d = D d + R d = I t o t d − I d en vervolgens de geplotte gegevens passen met de niet-lineaire functie

waar NS 0 en k zijn constanten die worden geschat door de niet-lineaire aanpassing. De functie wordt uitgedrukt in termen van alleen modelwaarden en wordt aangepast aan de curve van de gegevens. Dus, hebben verkregen NS van de niet-lineaire aanpassing, het aantal terugvorderingen R kan in de tijd worden beschreven door de simpele observatie dat het wordt gegeven door de geschaalde verhuizingen, Rm uit het SIR-model (1) minus het aantal doden, NS van Vergelijking (3),


Achtergrond

Deze serie is niet bedoeld om je snel een aantal plots met veel kleurrijke rondingen te laten zien die je daarvan moeten overtuigen mijn model kan coronavirusgevallen over de hele wereld perfect voorspellen. In plaats daarvan zal ik alle achtergronden uitleggen die nodig zijn om deze modellen te begrijpen, uw eigen mening over deze modellen te vormen en uw eigen ideeën te implementeren. Je hebt alleen wiskunde op middelbare schoolniveau nodig om de uitleg te volgen. Je hebt een goed begrip van python nodig om de programmeeronderdelen te volgen.

We willen infectieziekten modelleren. Deze ziekten kunnen zich van het ene lid van een populatie naar het andere verspreiden. We proberen inzicht te krijgen in hoe snel ze zich verspreiden, welk deel van een populatie ze besmetten, welk deel sterft, enz. Een van de gemakkelijkste manieren om ze te modelleren (en de manier waarop we zijn gericht op hier) is met een compartimentenmodel. Een compartimentenmodel verdeelt de populatie in verschillende compartimenten, bijvoorbeeld:

  • Svatbaar (kan nog steeds besmet zijn, “gezond”)
  • lbesmet
  • Rhersteld (waren al geïnfecteerd, kunnen niet opnieuw geïnfecteerd raken)

Dat wil zeggen, we kunnen een populatie hebben van N=1000 (bijvoorbeeld 1000 mensen) en we weten dat 400 mensen besmet zijn op tijdstip t (bijvoorbeeld t=7 dagen na het uitbreken van de ziekte). Dit wordt aangegeven met S(7) = 400. Het SIR-Model stelt ons in staat om, alleen door enkele initiële parameters in te voeren, alle waarden S(t), I(t), R(t) voor alle dagen t te krijgen. Ik zal nu de benodigde variabelen introduceren met een eenvoudig voorbeeld:

We hebben een nieuwe ziekte, ziekte X. Voor deze ziekte is de kans dat een besmet persoon een gezond persoon besmet 20%. Het gemiddelde aantal mensen waarmee een persoon per dag in contact is, is 5. Dus per dag ontmoet een geïnfecteerd persoon 5 mensen en infecteert ze elk met een kans van 20%. We verwachten dus dat deze persoon 1 persoon (20% ⋅ 5 = 1) per dag zal besmetten. Dit is β (“bèta”), het verwachte aantal mensen dat een besmette persoon per dag infecteert.

Nu kan men zien dat de aantal dagen dat een besmette persoon de ziekte heeft en kan verspreiden uiterst belangrijk is. We bellen dit nummer NS. Als D=7 loopt een besmette persoon zeven dagen rond om de ziekte te verspreiden en besmet hij 1 persoon per dag (omdat β=1). We verwachten dus dat een besmette persoon 1⋅7 (1 per dag maal 7 dagen) = 7 andere mensen infecteert. Dit is het basisreproductiegetal R₀, het totale aantal mensen dat een besmette persoon infecteert. We hebben zojuist een intuïtieve formule gebruikt: R₀ = β ⋅ D.

We hebben eigenlijk niets anders nodig, alleen een kleine notatie: γ ("gamma") zal 1/D zijn, dus als je denkt aan D als het aantal dagen dat een geïnfecteerde persoon de ziekte heeft, u kunt denken aan γ als de snelheid van herstel, of het percentage geïnfecteerde personen dat per dag herstelt. Als er nu bijvoorbeeld 30 mensen besmet zijn en D=3 (ze zijn dus drie dagen besmet), dan zal per dag 1/3 (dus 10) van hen herstellen, dus γ=1/3. Met γ = 1/D, dus D = 1/γ , en R₀ = β ⋅ D, volgt dat R₀ = β / γ.

Hier zie je de belangrijkste variabelen en hun definities nog eens terug:

  • N: totale populatie
  • NS): aantal mensen vatbaar op dag t
  • Het): aantal mensen besmet op dag t
  • R(t): aantal mensen hersteld op dag t
  • : verwacht aantal mensen dat een besmette persoon per dag infecteert
  • NS: aantal dagen dat een besmette persoon de ziekte heeft en kan verspreiden
  • : het aandeel geïnfecteerden dat per dag herstelt (γ = 1/D)
  • R₀: het totale aantal mensen dat een besmette persoon infecteert (R₀ = β / γ)

Hoe de veerconstante te berekenen

Er zijn twee eenvoudige benaderingen die u kunt gebruiken om de veerconstante te berekenen, met behulp van de wet van Hooke, naast enkele gegevens over de sterkte van de herstellende (of uitgeoefende) kracht en de verplaatsing van de veer vanuit zijn evenwichtspositie, of met behulp van de elastische potentiële energie vergelijking naast cijfers voor de arbeid die is verricht bij het verlengen van de veer en de verplaatsing van de veer.

Het gebruik van de wet van Hooke is de eenvoudigste manier om de waarde van de veerconstante te vinden, en je kunt de gegevens zelfs zelf verkrijgen via een eenvoudige opstelling waarbij je een bekende massa ophangt (met de kracht van het gewicht gegeven door ​F​ = ​mg​) van een veer en noteer de verlenging van de veer. Negeren van het minteken in de wet van Hooke (aangezien de richting er niet toe doet voor het berekenen van de waarde van de veerconstante) en delen door de verplaatsing, ​x​, geeft:

Het gebruik van de formule voor elastische potentiële energie is een even eenvoudig proces, maar het leent zich niet zo goed voor een eenvoudig experiment. Als u echter de elastische potentiële energie en de verplaatsing kent, kunt u deze berekenen met:

In ieder geval kom je uit op een waarde met eenheden van N/m.


Resultaten

De R0 verkregen voor de gegevens tot 19 april 2020, met behulp van geschatte SIR-modelparameters, was 1,02 [betrouwbaarheidsinterval (BI) van 0,75-1,29] met een RMSE van 7,72. De modelvoorspelling voor het totale aantal gevallen en het totale aantal werkelijk gerapporteerde gevallen wordt geïllustreerd in figuur 1.

Het gerapporteerde aantal totale gevallen (voor gegevens tot 19 april 2020) en het voorspelde aantal totale aantal gevallen door het model (RMSE = 7.72124)

De R0 berekend door het model met behulp van gegevens tot 30 april (RMSE van 172,44) was 1,66 (BI van 0,98-2,33). Figuur 2 illustreert de voorspelling van het totale aantal gevallen door het model in vergelijking met het totale aantal gemelde gevallen.

Het gerapporteerde aantal totale gevallen (voor gegevens tot 30 april 2020) en het voorspelde aantal totale aantal gevallen door het model (RMSE = 172,444)

De lage RMSE van het eerste model geeft aan dat het model meer representatief is voor de spreiding in Sri Lanka. De exponentiële groeisnelheidsmethode en de maximale waarschijnlijkheidsschattingsmethode gaven een R0 van 0,93 [betrouwbaarheidsinterval (BI) van 0,77-1,10] en een R0 van 1,23 (CI van 0,94-1,57), wanneer toegepast op deze dataset. Dit wordt geïllustreerd in Fig. 3.

R0-schattingen en 95%-betrouwbaarheidsintervallen van het SIR-model, exponentiële groeimethode (EG) en de maximale waarschijnlijkheids- (ML) schattingsmethode

Het real-time reproductiegetal (R) bleek een aanzienlijke variatie in de tijd te vertonen, variërend van 0,69 (95% geloofwaardige intervallen van 0,45–0,97) tot 2,20 (95% geloofwaardige intervallen van 1,65–2,83). Het is duidelijk dat implementatie van initiële beheersmaatregelen de overdraagbaarheid verminderde, die opnieuw steeg met de detectie van de twee clusters van gevallen. Ondanks deze tegenslag neemt de doorlaatbaarheid R echter weer af. Het dagelijkse totale aantal patiënten wordt geïllustreerd in Fig. 4 en 5 geven de variabiliteit van R in de tijd weer.

Het totale aantal dagelijks gemelde gevallen van 11 maart 2020 tot 7 mei 2020

De variatie van R, de overdraagbaarheid van COVID-19, van 11 maart 2020 tot 7 mei 2020


Modellering van infectieziekten

Het gebruik van wiskunde om de verspreiding van ziekten te modelleren is een ongelooflijk belangrijk onderdeel van de voorbereiding op mogelijke nieuwe uitbraken. As well as providing information to health workers about the levels of vaccination needed to protect a population, it also helps govern first response actions when new diseases potentially emerge on a large scale (for example, Bird flu, SARS and Ebola have all merited much study over the past few years).

The basic model is based on the SIR model – this is represented by the picture above (from Plus Maths which has an excellent and more detailed introduction to this topic). The SIR model looks at how much of the population is susceptible to infection, how many of these go on to become infectious, and how many of these go on to recover (and in what timeframe).

Another important parameter is R0, this is defined as how many people an infectious person will pass on their infection to in a totally susceptible population. Some of the R0values for different diseases are shown above. This shows how an airbourne infection like measles is very infectious – and how malaria is exceptionally hard to eradicate because infected people act almost like a viral storage bank for mosquitoes.

One simple bit of maths can predict what proportion of the population needs to be vaccinated to prevent the spread of viruses. The formula is:

waar Vt is the proportion of the population who require vaccinations. In the case of something like the HIV virus (with an R0 value of between 2 and 5), you would only need to vaccinate a maximum of 80% of the population. Measles however requires around 95% vaccinations. This method of protecting the population is called herd immunity


This graphic above shows how herd immunity works. In the first scenario no members of the population are immunised, and that leads to nearly all the population becoming ill – but in the third scenario, enough members of the population are immunised to act as buffers against the spread of the infection to non-immunised people.

The equations above represent the simplest SIR (susceptible, infectious, recovered) model – though it is still somewhat complicated!

dS/dt represents the rate of change of those who are susceptible to the illness with respect to time. dI/dt represents the rate of change of those who are infected with respect to time. dR/dt represents the rate of change of those who have recovered with respect to time.

For example, if dI/dt is high then the number of people becoming infected is rapidly increasing. When dI/dt is zero then there is no change in the numbers of people becoming infected (number of infections remain steady). When dI/dt is negative then the numbers of people becoming infected is decreasing.

The constants β and ν are chosen depending on the type of disease being modelled. β represents the contact rate – which is how likely someone will get the disease when in contact with someone who is ill. ν is the recovery rate which is how quickly people recover (and become immune.

ν can be calculated by the formula:

where D is the duration of infection.

β can then be calculated if we know R0 by the formula:

Modelling measles

So, for example, with measles we have an average infection of about a week, (so if we want to work in days, 7 = 1/ν and so ν = 1/7). If we then take R0 = 15 then:

Therefore our 3 equations for rates of change become:

Unfortunately these equations are very difficult to solve – but luckily we can use a computer program to plot what happens. We need to assign starting values for S, I and R – the numbers of people susceptible, infectious, recovered (immune) from measles. Let’s say we have a total population of 11 people – 10 who are susceptible, 1 who is infected and 0 who are immune. This gives the following outcome:

This shows that the infection spreads incredibly rapidly – by day 2, 8 people are infected. By day 10 most people are immune but the illness is still in the population, and by day 30 the entire population is immune and the infection has died out.

An illustration of just how rapidly measles can spread is provided by the graphic above. This time we start with a population of 1000 people and only 1 infected individual – but even now, within 5 days over 75% of the population are infected.

This last graph shows the power of herd immunity. This time there are 100 susceptible people, but 900 people are recovered (immune), and there is again one infectious person. This time the infection never takes off in the community – those who are already immune act as a buffer against infection.

If you enjoyed this post you might also like:

Differential Equations in Real Life – some other uses of differential equations in modelling predator-prey relationships between animal populations.

Essential resources for IB students:

Revision Village has been put together to help IB students with topic revision both for during the course and for the end of Year 12 school exams and Year 13 final exams. I would strongly recommend students use this as a resource during the course (not just for final revision in Y13!) There are specific resources for HL and SL students for both Analysis and Applications.

There is a comprehensive Questionbank takes you to a breakdown of each main subject area (e.g. Algebra, Calculus etc) and then provides a large bank of graded questions. What I like about this is that you are given a difficulty rating, as well as a mark scheme and also a worked video tutorial. Really useful!

The Practice Exams section takes you to a large number of ready made quizzes, exams and predicted papers. These all have worked solutions and allow you to focus on specific topics or start general revision. This also has some excellent challenging questions for those students aiming for 6s and 7s.

Each course also has a dedicated video tutorial section which provides 5-15 minute tutorial videos on every single syllabus part – handily sorted into topic categories.

I’ve put together four comprehensive pdf guides to help students prepare for their exploration coursework and Paper 3 investigations. The exploration guides talk through the marking criteria, common student mistakes, excellent ideas for explorations, technology advice, modeling methods and a variety of statistical techniques with detailed explanations. I’ve also made 17 full investigation questions which are also excellent starting points for explorations. The Exploration Guides can be downloaded here and the Paper 3 Questions can be downloaded here.


How can we calculate the constants of the SIR model? - Biologie

BME 332: Introduction to Biosolid Mechanics

Section 7: Constitutive Equations: Viscoelasticity

In the last section we introduced linear and nonlinear elastic constitutive models for materials and tissues. In many cases, elastic constitutive models work well when time dependent effects can be neglected. However, in those cases when time dependent effects cannot be neglected, we will need to utlize different constitutive models. Basically, time dependent effects indicate that the stress-strain behavior of a material will change with time. The classic material model for time dependent effects is viscoelasticity. As the name implies, viscoelasticity incorporates aspects of both fluid behavior (viscous) and solid behavior (elastic). In this section, we will cover physical characteristics for viscoelastic behavior as well as introduce basic mechanical analogs for describing viscoelastic behavior that can be used to derive constitutive equations. Finally, we will introduce general forms of a viscoelastic constitutive equation.

II. Characteristics of a Viscoelastic Material

Just like for elastic models, there are specific characteristics for viscoelastic models. These characteristics set viscoelastic models distinctly apart from elastic models. Most notably, we know that elastic materials store 100% of the energy due to deformation. However, viscoelastic materials do not store 100% of the energy under deformation, but actually lose or dissipate some of this energy. This dissipation is also known as hysteresis. Hysteresis explicitly requires that the loading portion of the stress strain curve must be higher than the unloading curve. Thus, from a stress-strain curve we can see the hysteresis as the area between the loading and unloading curve:

The ability to dissipate energy is one of the main reasons for using viscoelastic materials for any application to cushion shock, from running shoes to packing materials. The two other main characteristics associated with viscoelastic materials are stress relaxation and creep. Stress relaxation refers to the behavior of stress reaching a peak and then decreasing or relaxing over time under a fixed level of strain, as shown below:

Creep is in some sense the inverse of stress relaxation, and refers to the general characteristic of viscoelastic materials to undergo increased deformation under a constant stress, until an asymptotic level of strain is reached, as shown below:

Any materials that exhibit hysteresis, creep or stress relaxation can be considered viscoelastic materials. In comparison, elastic materials do not exhibit energy dissipation or hysteresis as their loading and unloading curve is the same. Indeed, the fact that all energy due to deformation is stored is a characteristic of elastic materials. Furthermore, under fixed stress elastic materials will reach a fixed strain and stay at that level. Under fixed strain, elastic materials will reach a fixed stress and stay at that level with no relaxation.

III. Mechanical Analogs for Viscoelastic Materials

The classic description and way to derive viscoelastic consitutive models is through the use of mechanical analogs. These are simple mechanical models for fluid and solid representations that are put together to produce viscoelastic effects. The simplest mechanical analog for a linear elastic material is a spring:

The simple constitutive relationship for a spring relates the force (and stress by extension when force is divided by area) to the elongation or displacement (and strain by extension when displacement is normalized by length of the spring):

where Fs is the spring force, E is the elastic modulus of the spring, and us represents the spring displacment. The mechanical analog for a Newtonian fluid is a dashpot. The simple constitutive relationship for a dashpot indicates that the force in the fluid depends on the rate the dashpot is displaced, or equivalently the velocity of the dashpot.

Also, the constitutive parameter that relates force (stress) to displacement rate (strain rate) is viscositeit:, which we denote as m . Thus, the constitutive equation for a fluid may be written as:

where the dot over the u in the equation indicates differentiation with respect to time and the superscript d denotes "dashpot".

By making various combinations of spring and dashpot models, we can simulate the behavior of a viscoelastic material, including stress relaxation and creep.

The simplest combination of the spring and dashpot is to put the spring in series with the dashpot:

This combination is known as the Maxwell model. As we will see, this model actually represents a fluid since it relaxes completely to zero stress and undergoes creep indefinitely. To derive the constitutive relation, we note first examine the kinematics of the model. It is clear from the geometry of the model that the total displacement will be the spring displacement plus the dashpot displacement:

However, the constitutive displacement for the dashpot is written in terms of the dashpot displacement differentiated with respect to time. Thus, to be consistent, we will differentiate the above expression for the total displacement with respect to time. This gives:

We need to determine both the displacement rate for the solid and the fluid. For the fluid, we simply rearrange the fluid constitutive equation:

To get the solid displacement rate, we first rearrange the constitutive equation for the spring:

we then differentiate this constitutive relationship with respect to time to obtain:

since E does not depend on time. Putting these together we get the total displacement rate as:

By normalizing the force by area and the displacement by length, we can get the analogous stress-strain relation as:

were the s and d superscripts have been dropped. We also have the following dimensions for the spring constant E and the viscosity m :

If we assume constants q1, p0 and p1, we can rewrite the above equation as:

where the constants are defined as:

Now, let us consider the behavior of a Maxwell material under conditions of both fixed strain and fixed stress. We assume that strain and stress can be applied instantaneously to a fixed level. In reality, instantaneous stress or strain cannot be reached, but we assume that in the limit, for example, if the test is run long enough, that the stress or strain can be considered to be reached instantaneously. We start with the response to an instantaneous strain that is then held constant over time:

In this case, the rate of change of strain is zero. Thus the equation for the Maxwell material becomes:

To solve for the stress under fixed strain for a Maxwell material, we need to solve the ordinary differential equation in time. To do this, we can again use the symbolic toolbox in MATLAB. First, we divide through by p1 to obtain:

We can then use the dsolve option in the MATLAB symbolic toolbox to solve the ordinary differential equation in time. The differential equation is put in parentheses after we define the symbols:

>> syms sig
>> syms sig p0 p1
>> dsolve('Dsig + (p0/p1)*sig = 0')

MATLAB returns a constant C1 which is an integration constant. We can write the equation as:

To solve for the integration constant C1, we need to determine the initial condition at time t = 0+ on the stress. To do this, we need to integrate the Maxwell constitutive equation from time t = 0- to some short time t = tn just after the strain is applied:

we know that e at t = 0- is 0 and e at t = tn is e 0. Also, s at t = 0- is 0 and s at t = tn is s 0. For the last term, we assume that tn is a small number. As we keep ramping up in shorter and shorter times, in the limit we approach tn = 0+, a nearly instantaneously applied strain. In dit geval,

This means the last term in the integration is zero and we are left with:

Thus, if we substitute in the initial condition to the solution of the Maxwell model constitutive equation we obtain the constant C1:

Thus, we have the general solution for the Maxwell model under a step strain as:

or if we write the constants in terms of the spring modulus and dashpot viscosity:

Now if we divide the stress as a function of time by the initial strain, we obtain the stress relaxation function G(t) for the Maxwell model:

This is a characteristic function that tells us how the stress relaxes for a given material. If we plot the the stress relaxation versus time we obtain the general curve:

Note that the stress completely relaxes out over time. This is actually characteristic of a viscoelastic fluid rather than a viscoelastic solid. Finally, we note that the ratio m /E gives us a dimension of time. This is a characteristic time of the material denoted as the relaxation time, and defined as:

Using the relaxation time, we can rewrite the relaxation function as:

Now let us consider the response of the Maxwell model to a unit step stress:

In this case, a stress s 0 is applied "instantaneously" and then kept constant over time. Thus, the rate of change of stress is zero and the Maxwell constitutive model becomes:

We can solve the above equation simply by integrating with respect to time. We can also do this in MATLAB as:

>> syms eps sig0 p0 q1
>> dsolve('Deps = (p0/q1)*sig0')

where again C1 is a constant of integration that we need to determine from initial conditions. Given the initial conditions derived before, we obtain:

This gives the strain versus time under a unit step stress as:

Similar to the stress relaxation, we can define a creep function J(t) by dividing the strain versus time reponse by the initial unit step stress. This gives the creep function Jas:

Note that the above equation will give a constantly increasing strain as a function of t:

This is also not characteristic of a viscoelastic solid, since we expect that the creep will reach an asymptotic level after a certain time.

Although the Maxwell model did not give us a realistic result for a viscoelastic solid (it is more representative of a viscoelastic fluid) it did allow us to illustrate some important aspects of a viscoelastic material including the stress relaxation function G(t) and the creep function J(t). We next consider another mechanical analog for a viscoelastic material, the Voight model, also known as the Kelvin-Voight model. Its geometry is shown below:

Again, we can make observations based on the geometry of the model. First, we note that the dashpot will constrain the spring to have the same deformation:

Second, we note that the total force F in the Voight model will be equal to the force in the dashpot plus the force in the spring:

By analogy, we can write a stress-strain differential equation as:

The above equation again illustrates an important characteristic of viscoelastic materials, namely that the stress in the material depends not only on the strain, but also on the strain rate.

Let's now see how the Voight model responds to a unit step stress and strain. First, we consider the unit step strain again:

For the last term in the integral, we have:

For the second term, we have:

If we look at the first term, it will give us the area under the stress versus time curve. This area must have a limit, otherwise we would have infinite stress. Thus, in the limit as the time becomes very short, the area must have a finite value. This is given by the dirac delta function (not to be confused with the kronecker delta), which has the property:

Thus, we have the following initial condition for a step strain test:

which gives us the following stress versus time response:

This gives a stress relaxation function as:

This function gives instantaneous stress relaxation due to the presence of the dashpot.

Next, let us consider the response of the Voight model to a unit step stress:

We again solve the differential equation in time for the Voight model using MATLAB. To get the equation into a form we can use for MATLAB, we first divide through by the viscosity and recognize that the stress is constant s 0 to get:

We can then solve this equation in MATLAB:

>> syms sig0 eps emod mvis
>> dsolve('Deps + (emod/mvis)*eps = sig0/mvis')

where again C1 is a constant of integration. To determine the constant C1, we again need to use the initial condition. In this case, under an instantaneously applied stress, the dashpot cannot move instantaneously. Thus, this means that there will be no displacement and therefore the strain at time t = 0 will be zero. Pluggin this initial condition into the above equation gives:

This give us the solution:

If we divide the above by s 0, we obtain the creep function J:

The above function gives us a strain versus time plot as shown below:

Note that no instantaneous elastic deformation is possible due to the restraint of the dashpot.

Note that if we rearrange the above equation we can write the strain as a function of time and the initial stress as:

Note that as mentioned before the dashpot does not allow instantaneous deformation to occur, but overtime the displacement creeps to an asymptotic level.

NS. Generalized Viscoelastic Constitutive Models

We can see from the mechanical analogs investigated in the previous section that the strain behavior over time of a viscoelastic material is a function of the creep function and the stress, while the stress behavior over time is a function of the stress relaxation function and the strain. Boltzmann (1844-1906) first generalized these observations by saying that for a simple bar subject to a stress s (t), that the increment in stress over a small time interval d t would be:

This assumes that the stress is continuous and differentiable in time. Given that the stress is related to the stress via the creep function, Boltzmann postulated that an increament of strain d e , which depends on the complete stress history up to time t, would be related to the increment of stress d s at the specific time increment from t to t through the creep function J at the time t - t as:

The complete strain at a time t would then be obtained by integrating the strain increments from time 0 to time t, over all the increments d t :

We can also make the same argument for an increment of stress d s through time t depending on the increment of strain d e at time t and the stress relaxation over the time t - t as:

Again, if we integrate over the entire strain history we obtain:

The above constitutive relationships can be generalized to three dimensions in tensorial form for a strain history from time t = - infinite to t as:

where now there is a tensorial stress relaxation function Gijkl.

We can write a similar relationship for strain in terms of a prescribed stress history and a tensorial creep function Jijkl as:

If we know that the loading starts at t = 0 and the stress and strain are zero at this point, we may write the above constitutive relationships as:

V. Fung's quasi-linear viscoelasticity theory

The constitutive equations presented in the mechanical analog section and the generalization seen above work for materials undergoing small strain. However, we have seen that most soft tissues undergo large deformation. To account for this, Fung (see Y.C. Fung, Biomechanics: Mechanical Properties of Living Tissues) introduced the concept of quasi-linear viscoelasticity (QLV) theory. In this case, Fung introduced a separation of the stress relaxation function into a time dependent part and an elastic portion. This gives:

where now the total stress relaxation function K depends on both time and the stretch ratio l in terms of a reduced relaxation function G that only depends only on time and an elastic response Te that depends on the stretch ratio. In this way, the theory is linear in the relaxation response but still accounts for nonlinear large deformation.

The stress response for an infinitesimal change in stretch can then be written as:

where T is the 1st PK stress as a function of time, and Te is an energy function. If we assume that superposition holds for the linear case, we can integrate the above equation to obtain:

Note that the form above is very similar to the linear viscoelastic model for stress, with Te the elastic response depending on the stretch ratio replacing the small strain rate.

Fung also proposed a general stress relaxation function for use in QLV. This reduced relaxation function is given as:

where c is a constant, t 1 and t 2 are characteristic times during relaxation, and E1 is the exponential integral:

note that as t goes to inifinity the above integral goes to zero and the reduced relaxation function reduces to:

QLV is an approximation used quite often for viscoelastic constitutive models of soft tissues.

VI. Example of Computing Stress Versus Time for a Quasilinear Viscoelastic Model

To give an example of a quasilinear viscoelastic model, consider the model proposed by Toms et al. (2002) in the Journal of Biomechanics for periodontal ligament. The periodontal ligament is a thin piece of soft tissue interposed between a tooth and the surrounding alveolar bone. The PDL provides cushioning between the tooth and the bone, as well as being a conduit for nutrients. It is also known that the mechanical behavior of the PDL is nonlinear and viscoelastic. To characterise this behavior, Toms et al. chose to use the quasilinear viscoelastic approach of Fung, where the general constitutive equation is represented by:

although they did not use the reduced relaxation function proposed by Fung. Instead they chose a general decaying exponential function of the form:

where e is the exponential function, t is time, and a,b,c,d,g,h are all coefficients to be determined by experiment. This was a function that had been previously proposed by Will et al. in 1972 for characterising the viscoelastic behavior of the PDL. For the instantaneous nonlinear elastic response, Toms et al chose a very commonly used form:

where T is the 1st PK stress, e is the exponential function, l is the principal stretch ratio, and A and B are constants to be determined experimentally. The strain energy function that would correspond to this instantaneous nonlinear elastic response can be derived by integrating the 1st PK stress over l . We can do this quite easily using the MATLAB symbolic toolbox as follows:

As always, we first define the symbols we will be using:

We then perform the integration in MATLAB using the int command, where we specify the function followed by l as the variable of integration:

which gives the strain energy function defined as:

Thus, the elastic response becomes the second derivative of the strain energy function with respect to the principal stretch ratio which gives us the 1st PK stress.

Let's consider a ramp displacement of the PDL with the following time constant:

Now let us calculate the 1st PK stress at some time past t0. To do this, Toms et al used the following test set up to determine the nonlinear elastic response as well as the time dependent stress decay:

which was used due to the thinness of the PDL and the difficulty in separately the PDL from the tooth and the alveolar bone. Based on Toms results, let's calculate a typical stress time reponse analytically. To do this, we write the general form of the quasilinear viscoelastic function:

where E is the large strain tensor. For the above ramp strain, we know that the strain rate after t = t0 is zero since the strain is constant. Therefore, we need only integrate to t0. Furthermore, we know that the strain rate up to t0 is a constant g , which we can substitute into the above equation to obtain:

we can now integrate a general expression for the 1st PK stress analytically in MATLAB, in similar fashion as to how we integrated the stress to get the strain energy function. We first define the symbols:

>>syms a b c d g h A B lam t tau t0 gam

we next write the function and integrate the function over t from 0 to t0:

MATLAB gives us the answer:

Now let us use the above analytical expression to compute stress versus time. We first enter values for all the constants in MATLAB: