Een semi-stochastisch cell-based model voor de ontwikkeling van een longcarcinoom

Technische Universiteit Delft Faculteit Elektrotechniek, Wiskunde en Informatica Delft Institute of Applied Mathematics

Een semi-stochastisch cell-based model voor de ontwikkeling van een longcarcinoom

Verslag ten behoeve van het Delft Institute of Applied Mathematics als onderdeel ter verkrijging

van de graad van

BACHELOR OF SCIENCE in TECHNISCHE WISKUNDE

door

RUNE VAN DER MEIJDEN Delft, Nederland Februari 2013

Copyright

© 2013 door de TU Delft. Alle rechten voorbehouden.

BSc verslag TECHNISCHE WISKUNDE

“Een semi-stochastisch cell-based model voor de ontwikkeling van een longcarcinoom”

RUNE VAN DER MEIJDEN

Technische Universiteit Delft

Begeleider Dr.ir. F. J. Vermolen

Overige commissieleden Dr.ing. D. Jeltsema

Drs. E.M. van Elderen

Dr. J.L.A. Dubbeldam

Februari, 2013

Delft

“If people do not believe that mathematics is simple, it is only because they do not realize how complicated life is” - John Louis von Neumann

Voorwoord “Ik woon op huisnummer 23 en ik werk in gebouw 23. Mijn geboortedatum is 23-03-1991. 23 is 0 + 3 + 1 + 9 + 9 + 1. Mijn naam bestaat uit de 18e, 21e, 14e en 5e letter van het alfabet. En 1418−5 + 21 is ook 23. Mijn bankrekeningnummer is 14.28.39.493. Uiteraard is (1-4) + (2 x 8) + (3 + 9) + (4 - 9 + 3) gelijk aan 23. En van mijn pincode kun je ook 23 maken. Denk ik.” Al sinds ik mij kan heugen ben ik gefascineerd door getallen en het feit dat overal een structuur in te vinden is. In de tijd dat ik fanatiek wielrenner was, kon ik bij het ontstaan van een kopgroep, aan de hand van de rugnummers van mijn medevluchters al snel concluderen of we de finish zouden halen of voortijdig teruggepakt zouden worden. Hiervoor hoefde ik alleen maar een gemeenschappelijke deler te vinden of te concluderen dat alle rugnummers priemgetallen waren. Als ik door een combinatie van optellen, aftrekken, vermenigvuldiging en deling van alle andere rugnummers het mijne kon maken stond ook de winnaar al vast. Uiteraard geloofde ik zelf niet ´echt in deze “meta-mathematische” samenhang tussen rugnummers en de uitslag van een wedstrijd, maar het leidde me wel af van de pijn in mijn benen en liet de 200km durende wedstrijden sneller voorbij gaan. In die zin zorgde het indirect dus wel voor betere prestaties. Wiskunde was mijn doping. Deze fascinatie voor getallen en structuur leidde er toe dat ik er voor koos mij toe te wijden aan de meest objectieve wetenschap der wetenschappen, de kunst van ‘het zeker weten’; ik besloot de bacheloropleiding Technische Wiskunde te gaan doen aan de TU Delft, een keuze waar ik tot op de dag van vandaag geen moment spijt van heb gehad. In de drie jaar die mijn bachelor heeft geduurd, heb ik een prachtige kijk in de wereld van zowel de abstracte als de concrete wiskunde gekregen. Dit project is het eindproduct van deze jaren. De kunst van het puzzelen met de meest fundamentele bewijzen is een prachtige uitdaging, maar nog meer voldoening haal ik uit het ontrafelen van complexe en maatschappelijk relevante problemen, om aan de hand van modellen en simulaties een antwoord te geven op vragen die aanvankelijk onbeantwoordbaar leken. Dit is de reden geweest mijn Bachelor eindproject te doen bij de afdeling Numerieke Methoden, bij begeleider Dr.Ir. Fred Vermolen. Het veld waarin hij werkzaam is - de biologische mathematica - is een onderwerp dat mij enorm fascineert. Ik wil hem bedanken voor zijn ondersteuning en de vrijheid die hij mij heeft geboden mijn eigen invulling aan het project te geven. Daarnaast wil ik ook zijn collega Dr. Amit Gefen van de universiteit van Tel Aviv bedanken voor zijn kennis van de meer cel-biologische aspecten. Tot slot wil ik u veel plezier wensen in het lezen en beoordelen van dit verslag.

Samenvatting In de geneeskunde wordt het gebruik van wiskundige modellen en computersimulaties steeds belangrijker. Zo kan steeds beter voorspeld worden hoe complexe systemen in het lichaam reageren op externe factoren zoals medicijnen. Een toepassing hiervan vindt men in kankeronderzoek. Cellen in het menselijk lichaam kunnen op twee manieren worden aangezet tot beweging. Actieve migratie wordt veroorzaakt doordat cellen elkaars aanwezigheid voelen en naar elkaar toe bewegen. Passieve migratie ontstaat als twee cellen tegen elkaar aan komen te liggen en door een afstotende kracht uit elkaar gaan bewegen. De fysieke eigenschappen van de cel en zijn omgeving be¨ınvloeden de mate waarin deze vormen van migratie plaatsvinden. Cellen maken ook een cyclus van groei en deling door, waarin een passieve “pauzestand” aan kan worden genomen als weefselgroei niet gewenst is. Als er een teveel aan cellen is, kunnen individuele cellen overgaan tot geprogrammeerde sterfte, de zogenaamde apoptose. De kansen dat een cel overgaat tot deling of sterfte worden be¨ınvloed door de mate waarin aanliggende cellen druk uitoefenen. Bij tumorcellen is dit mechanisme, dat ervoor zorgt dat cellen bij het juiste weefselvolume in “pauzestand” raken, ontregeld en gaan ze als gevolg hiervan ongecontroleerd en snel delen. Hierdoor wordt gezond weefsel weggedrukt en kan dit afsterven. Als een tumor zich uitbreidt naar ander weefsel spreken we van uitzaaiingen en lijdt het betreffende organisme aan kanker. In dit project is de manier waarop cellen elkaar onderling aanzetten tot migratie, groei en sterfte gemodelleerd en zijn de beginstadia van de ontwikkeling van een kolonie tumorcellen ge¨ımplementeerd. In het bijzonder kijken we hierbij naar ´e´enlagig plaveiselepitheel, zoals dat zich in de longblaasjes bevindt. Dit is gemodelleerd als een celkweek van ´e´en laag dik dat zich bevindt op een substraat dat ingekapseld is door een weefselwand. De beweging en groei worden beide met een Euler voorwaartse methode ge¨ıntegreerd. De nadruk is gelegd op een realistische en tijdseffici¨entie implementatie. Dit heeft geleid tot een model dat klaar is voor verder gebruik in onderzoek naar de eerste stadia van de ontwikkeling van een tumor.

Verklarende woordenlijst Begrip Epitheel Epitheelcellen Celmigratie Celcyclus Mitose Apoptose Substraat Stijfheid Vervormbaarheid

Omschrijving Weefsel dat het vrije oppervlak van het lichaam bekleedt Cellen die deel uit maken van het epitheel De manier waarop cellen binnen een weefsel bewegen onder invloed van elkaar De levenscyclus van cellen Fase van de celcyclus waarin deling plaatsvindt Geprogrammeerde celsterfte Voedingsbodem waarop cellen kunnen worden gekweekt en zich voortbewegen Mate waarin weefsel een tegenkracht geeft bij vervorming Mate waarin weefsel vervormd bij uitoefenen van kracht

Lijst van gebruikte symbolen Symbool Λ Ω α β γ δ  λ µ σ τ2 φ ψ E F M R T V cf h r b u b z

Omschrijving Celdruk Leefgebied Bewegingsparameter Mobiliteitsparameter Groeiparameter Mate van indeuking dochter cellen na celdeling Minimaal detecteerbare rek-energiedichtheid Uitdovingsco¨efficient Groeisnelheid van cel-straal Groeivariantie Verdubbelingstijd Fase van celcyclus Mitose-parameter Elasticiteitsmodulus Trekkracht van een cel op de ondergrond Rek-energiedichtheid Radius Duur van celfase Absolute bewegingssnelheid van een cel Wrijvingsco¨efficient Indeuking Locatie Richting van celdeling Bewegingsrichting van een cel

Eenheid (indien van toepassing) kP a m3 Ns s−1

m/s s Pa N Pa m s m/s m m m m/s

Inhoud Voorwoord

i

Samenvatting

iii

Verklarende woordenlijst

v

Lijst van gebruikte symbolen 1 Introductie 1.1 Belang van het onderzoek 1.2 Probleemstelling . . . . . 1.3 Opbouw verslag . . . . . . 1.4 Gebruikte software . . . .

vii

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

2 Inleiding in de Celbiologie 2.1 De cellen van het lichaam . . . . . . . 2.2 Epitheel . . . . . . . . . . . . . . . . . 2.3 Celmigratie . . . . . . . . . . . . . . . 2.3.1 Botsende cellen . . . . . . . . . 2.3.2 Communicatie tussen cellen via 2.4 Het leven van een cel . . . . . . . . . . 2.4.1 Celcyclus . . . . . . . . . . . . 2.4.2 Apoptose . . . . . . . . . . . . 2.5 Tumorvorming en kanker . . . . . . . 2.5.1 Oorzaak . . . . . . . . . . . . . 2.5.2 Longkanker . . . . . . . . . . . 2.5.3 De tumorcel . . . . . . . . . . . 2.5.4 Tumorgrootte . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

1 1 1 2 2

. . . . . . . . . . . . . . . . . . . . . . . . substraat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

3 3 3 4 4 5 6 6 7 7 8 9 9 10

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

11 11 12 13 13 15 16 17 17 18 18 20

. . . .

3 Mathematisch model 3.1 Epitheelweefsel in de longen . . . . . . . 3.2 Beweging van cellen . . . . . . . . . . . 3.2.1 Bewegingsvergelijking . . . . . . 3.2.2 Actieve migratie . . . . . . . . . 3.2.3 Passieve migratie . . . . . . . . . 3.2.4 Mitotische cellen . . . . . . . . . 3.2.5 Bewegingsrichting . . . . . . . . 3.2.6 Beweging numeriek ge¨ıntegreerd 3.3 Cel cyclus . . . . . . . . . . . . . . . . . 3.3.1 Groei . . . . . . . . . . . . . . . 3.3.2 Checkpoint . . . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . . . . . . . . .

Inhoud . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . van tumorcellen . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

21 22 23 23 23 23 24 25 25 25

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

27 27 27 28 28 28 28 29 29 30 30 30

5 Resultaten 5.1 Migratie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Starttoestand eindsimulatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Simulatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32 32 33 34

6 Conclusie

36

7 Discussie 7.1 Aanbevelingen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37 38

3.4

3.5

3.3.3 Mitose . . . . . . . . . . . . 3.3.4 Apoptose . . . . . . . . . . Invloeden van externe druk . . . . 3.4.1 Celdruk . . . . . . . . . . . 3.4.2 Celcyclus . . . . . . . . . . 3.4.3 Drukafhankelijke sterfte . . Tumorvorming . . . . . . . . . . . 3.5.1 De tumorcel . . . . . . . . . 3.5.2 Intercellulaire communicatie 3.5.3 Celcyclus van een tumorcel

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . en migratie . . . . . . .

4 Numerieke implementatie 4.1 Gebruikte Numerieke methode . . . . 4.1.1 Beweging . . . . . . . . . . . . 4.1.2 Deling . . . . . . . . . . . . . . 4.1.3 Celgroei . . . . . . . . . . . . . 4.2 Implementatie . . . . . . . . . . . . . . 4.2.1 Vectorisatie . . . . . . . . . . . 4.2.2 Celcyclus . . . . . . . . . . . . 4.2.3 Sterfte en deling . . . . . . . . 4.3 Parameterkeuze . . . . . . . . . . . . . 4.3.1 Morfologie en fysica van de cel 4.3.2 Levenscyclus van de cel . . . .

Literatuurlijst

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

I

A Appendix III A.1 Implementatie in matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III A.1.1 Simulatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III A.1.2 Subfuncties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VIII

ix

Hoofdstuk 1

Introductie In dit project is onderzoek gedaan naar de mogelijkheden om de beginstadia van de ontwikkeling van een tumor realistisch te modelleren en effici¨ent te implementeren. Hieronder zal kort de probleemstelling worden toegelicht en de opbouw van het verslag worden uitgelegd. Ook zal de gebruikte software aan de orde komen.

1.1

Belang van het onderzoek

Zoals eerder vermeld, wint het gebruik van geavanceerde wiskundige modellen terrein in de medische wetenschap. Dit komt door een combinatie van betere modellen, snellere rekenmethoden en toenemende rekenkracht. Het gebruik van deze modellen zorgt ervoor dat vooraf voorspellingen kunnen worden gedaan over de uitkomst van een experiment of de werking van een medicijn. Hierdoor kan het aantal medische experimenten drastisch worden teruggedrongen. Dit spaart kosten, versnelt de ontwikkeling van nieuwe medicijnen en vermindert het belang van dierproeven. Het is hiervoor wel belangrijk dat de modellen goed en de simulaties tijds-effici¨ent zijn. Voor dit specifieke onderzoek is het interessant te zien hoe een weefsel op cel-niveau beweegt en hoe deze dynamiek verandert als een tumorcel uitgroeit tot een kolonie. Hierdoor kan meer inzicht worden vergaard in de manier waarop de eerste stadia van de ontwikkeling van een tumor zich manifesteren. Dit kan van nut zijn voor onderzoek op zowel celbiologisch als geneeskundig vlak.

1.2

Probleemstelling

In dit project is de probleemstelling aanvankelijk redelijk open gebleven. Het is een proces geweest waarbij hernieuwde kennis vanuit literatuur zorgde voor iteraties in zowel het model als de implementatie hiervan. Hierdoor is uiteindelijk gekozen om de focus van dit project te houden op het zo realistisch mogelijk maken van het model en het zo tijdseffici¨ent mogelijk implementeren hiervan. Hierdoor is de hoofdvraag die in dit onderzoek beantwoord is: “Is de ontwikkeling van een tumorkolonie realistisch te modelleren op celniveau?” Belangrijke deelvragen hierbij zijn: ˆ Hoe vindt intercellulaire communicatie plaats en hoe regelt dit de migratie van cellen? ˆ Wat is de invloed van de druk van omliggende cellen op de celcyclus, inclusief groei, deling en sterfte?

Introductie ˆ Wat is het verschil in fysiologische eigenschappen tussen gezonde cellen en tumorcellen? ˆ Hoe kunnen eigenschappen als verdubbelingstijd van een type weefsel en de duur van de celcyclus meegenomen worden in dit model? ˆ Wat is de juiste methode om dit model tijdseffici¨ent en toch stabiel en nauwkeurig te implementeren?

1.3

Opbouw verslag

Het verslag is opgebouwd uit een viertal hoofdstukken en wordt afgesloten door een conclusie en een discussie met aanbevelingen. Hoofdstuk 2 behandelt de celbiologische achtergrond die nodig is voor het opstellen van het model. Hierin wordt beschreven hoe cellen migreren onder invloed van elkaar, hoe hun celcyclus verloopt, hoe een tumor ontstaat en hoe de beweging en levensloop van tumorcellen zich verhoudt tot die van gezonde cellen. In hoofdstuk 3 wordt allereerst kort het oorspronkelijke model voor cel-migratie, zoals beschreven in [1], behandeld en wordt vervolgens dieper ingegaan op de uitbreiding hierop voor tumorvorming. Hierbij komen zowel de celcyclus als de implementatie van de tumor aan de orde. Ook zullen de keuze van delings- en sterfteparameters op basis van kansrekening worden afgeleid. In hoofdstuk 4 komt vervolgens de numerieke implementatie van dit model aan de orde. Eerst wordt hierin de gebruikte numerieke methode toegelicht. Daarna volgt de implementatie van beweging, deling en sterfte en wordt de vectorisatie behandeld. Tot slot wordt de keuze van de gebruikte parameters gegeven. De Matlab-implementatie is te vinden in Appendix A. In hoofdstuk 5 worden de belangrijkste resultaten weergegeven en beschreven. In hoofdstuk 6 worden de hierbij horende conclusies getrokken en in hoofdstuk 7 volgt een discussie hierop en worden aanbevelingen voor verder onderzoek gedaan.

1.4

Gebruikte software

In dit project zijn twee softwarepakketten gebruikt. Beide zijn standaard binnen de bachelor en zullen hieronder daarom slechts zeer beknopt worden toegelicht. ˆ Alle simulaties in dit project zijn ge¨ımplementeerd in de technische softwareomgevind Matlab ©(R2012b, The MathWorks, Natick, USA). De gebruikte .m-files zijn bijgevoegd in de bijlage. ˆ Alle tekstverwerking is gedaan in LATEX, het meest gebruikte pakket macro’s voor opmaaktaal TEX.

2

Hoofdstuk 2

Inleiding in de Celbiologie In dit project wordt een kolonie van dierlijke cellen gemodelleerd en gesimuleerd hoe hier een tumor in groeit. In het kader van onderzoek naar longkanker, wordt hier gekeken naar epitheelcellen die zich in de longblaasjes bevinden. Deze cellen hebben bepaalde fysieke eigenschappen zoals vorm en grootte. Ook is er interactie tussen deze cellen onderling. Ze bewegen naar elkaar toe doordat ze elkaars aanwezigheid kunnen voelen en als cellen botsen stoten ze elkaar weer af. Deze beweging van cellen wordt de celmigratie genoemd. Iedere cel maakt een levenscyclus - de celcyclus- door, waarin ze groeien en zich vervolgens kunnen splitsen in twee nieuwe cellen. Als er te veel cellen zijn, of als cellen niet meer nodig zijn en in de weg liggen, kunnen deze vernietigd worden. Deze geprogrammeerde celsterfte heet apoptose. Als het DNA van een cel muteert kan het zijn dat er in dit proces van celdeling en -sterfte iets mis gaat en kunnen de cellen ongecontroleerd snel gaan groeien en delen, wat leidt tot vorming van een tumor. Het betreffende organisme lijdt dan aan kanker. In het geval van de in dit project gesimuleerde situatie spreken we van longkanker. Er volgt nu een meer inhoudelijke introductie in de relevante celbiologie en de hierboven genoemde begrippen zodat duidelijk wordt wat er gesimuleerd wordt, waar rekening mee moet worden gehouden en waarin het uiteindelijke model overeenkomt met en verschilt van de werkelijkheid.

2.1

De cellen van het lichaam

Het menselijk lichaam bestaat uit een groot aantal cellen. Schattingen voor een gemiddeld persoon lopen uiteen, maar bevinden zich in de orde van 10-100 triljoen (1013 −1014 ) individuele cellen. Over het aantal verschillende typen cellen dat in het menselijk lichaam voorkomt is enige discussie. Afhankelijk van de definitie, bevindt dit getal zich tussen de 210 [2] en 411 [3]. Deze cellen verschillen van elkaar in fysieke eigenschappen en in functie. Voorbeelden hiervan zijn rode bloedcellen die zuurstof transporteren en zenuwcellen die signalen van en naar de hersenen geleiden. Een ander voorbeeld zijn de epitheelcellen. Deze cellen bevinden zich op de oppervlakten van organen, bijvoorbeeld de longen, nieren en huid.

2.2

Epitheel

De oppervlakten van het lichaam en holtes in het lichaam worden bedekt met epitheelweefsel. Dit weefsel bestaat uit aan elkaar grenzende cellen van ´e´en of meerdere lagen dik. Tussen deze cellen kunnen tight junctions worden gevormd, verbindingen die de cellen bij elkaar houden. De manier waarop deze cellen tegen elkaar aan liggen zorgt er voor dat bepaalde stoffen wel, en andere stoffen niet worden doorgelaten. Hierdoor kan het weefsel verschil-

Inleiding in de Celbiologie lende functies tegelijk hebben. Zo bestaat de menselijke huid uit een dikke laag epitheelcellen die een barri`ere vormt om enerzijds het inwendige van het organisme te beschermen tegen schadelijke externe factoren. Anderzijds laat de huid vocht en afvalstoffen door en kan hij voedingsstoffen opnemen. Ondanks dat cellen uit epitheelweefsel veel kunnen verschillen van vorm en de weefsels onderling sterk vari¨eren in dikte, zijn er wel een aantal gemeenschappelijke eigenschappen. Het weefsel is de afscheiding van het lichaam met de buitenwereld. Het grenst dan ook met ´e´en kant aan het uitwendige, zoals bijvoorbeeld de lucht in de longen of het voedsel in de darmen. De andere kant grenst aan het inwendige, meestal aan bindweefsel van bijvoorbeeld een bloedvat. Het wordt gescheiden van dit weefsel door een laagje extracellulaire matrix, genaamd de basale lamina.

Figuur 2.1: Typen epitheelcellen. ‘Simple squamous’ is eenlagig plaveisel. [4] In dit project is gekeken naar het celweefsel van een longblaasje. Dit is een eenlagig plaveiselcelepitheel (zie figuur 2.1), dat een hoge doorlatendheid voor gas kent. Zuurstof kan hierdoor makkelijk worden opgenomen door de longen en CO2 kan even eenvoudig weer worden uitgestoten. We bekijken een kolonie van cellen die in een laboratorium is gekweekt. De cellen bevinden zich op een substraat, de voedingsbodem waarop de cellen zijn gekweekt en waarop ze zich voort kunnen bewegen. Deze is te vergelijken met de basale lamina van het (long)epitheel.

2.3

Celmigratie

Een belangrijk fenomeen binnen de celbiologie is celmigratie, waarbij cellen zich actief verplaatsen door fysische of chemische prikkels van buitenaf, zoals beschreven staat in [7] en weergegeven in figuur 2.2. Deze signalen worden opgevangen door bepaalde eiwitten in het celmembraan. Deze voelen gezamenlijk de gradi¨ent van de signaalsterkte, waarna de cel zijn looprichting bepaalt. Dit kan naar de gradi¨ent toe, of juist er van af zijn, afhankelijk van of het signaal een gunstige of ongunstige omstandigheid voor de cel betekent.

2.3.1

Botsende cellen

Als twee cellen naar elkaar toe bewegen is de stijfheid van de ondergrond waar de cellen op bewegen bepalend voor wat er gebeurt als de cellen tegen elkaar aan botsen. De stijfheid 4

Figuur 2.2: Celmigratie. [8]

Celmigratie is de mate waarin het substraat een tegenkracht geeft als hij wordt ingedeukt. Een hoge stijfheid, betekent een hoge tegenkracht. De elasticiteitsmodulus van Young (E) is een mate om deze stijfheid te kwantificeren, de eenheid hiervan is Pa. Uit laboratoriumonderzoek met endotheelcellen [9] blijkt dat voor een ondergrond met een lage elasticiteitsmodulus (E < 2500 Pa) de cellen tegen elkaar aanbotsen en vervolgens tegen elkaar aan blijven liggen. Voor een grote modulus (E > 5500 Pa) geldt dat de cellen tegen elkaar aan botsen en vervolgens weer van elkaar af lopen. Voor tussenliggende waarden blijven de cellen tegen elkaar aan ‘botsen’: een herhaald naar elkaar toe en van elkaar af migreren, waarbij ze dicht bij elkaar blijven, mits de omstandigheden niet veranderen.

2.3.2

Communicatie tussen cellen via substraat

De observatie dat cellen op een middelmatig stijf substraat van elkaar en naar elkaar toe blijven bewegen, heeft tot de hypothese geleid dat cellen via het substraat kunnen communiceren. In onderzoek [9] is dit aangetoond met behulp van trekkracht-microscopie.

Levende cellen voeren een trekkracht uit op het substraat. De dit leidt tot een spanning (trekkracht per oppervlakte) die het substraat vervormd naar buiten toe. Dit mechanische signaal kan opgepikt worden door nabijgelegen cellen die als reactie hierop richting de ‘trekkende’ cel bewegen. Cellen oefenen een hogere trekkracht uit op een stijf substraat dan op een zwakker substraat. Echter vervormt een zwak substraat makkelijker, waardoor de totale vervorming toch groter is en minder snel uitdooft van de bron af. Omliggende cellen zullen dit signaal dus vanaf grotere afstand oppikken. Figuur 2.3: Substraatdeformatie. [10] Dichtbij elkaar gelegen cellen voelen dus elkaars aanwezigheid en beginnen naar elkaar toe te ‘lopen’ tot ze tegen elkaar aan botsen en vervolgens - afhankelijk van de stijfheid van het substraat- tegen elkaar aan blijven liggen, tegen elkaar aan blijven stuiteren of van elkaar af bewegen. In een kolonie van cellen is dit een complex systeem, waarin cellen de trekkracht van een hoop cellen tegelijk voelen. Deze trekkracht kan beschouwd worden als een potentiaal, en dus kunnen deze krachten bij elkaar opgeteld worden. De richting waarin de cel loopt wordt bepaald door alle omliggende cellen samen. De snelheid waarmee de cel beweegt zal nagenoeg constant blijken. De formules waarmee dit proces van vervorming van het substraat, het waarnemen hiervan en de mate waarin dit de looprichting van de cel bepaalt en wat er gebeurt als twee cellen tegen elkaar aan botsen, komen Hoofdstuk 3 aan bod. 5

Inleiding in de Celbiologie

2.4

Het leven van een cel

Alle eukaryotische cellen - cellen met een kern - maken biologisch gezien dezelfde reeks van levensfasen door die uiteindelijk leidt tot de deling van een moedercel in twee dochtercellen. Dit proces heet de celcyclus. Het delen van ´e´en cel in twee volwaardige dochtercellen heet mitose. De verschillende fasen van de mitose zullen kort toegelicht worden, zodat later gemaakte keuzes in het model verdedigd kunnen worden. Meiose, wat alleen plaatsvindt bij mannelijke organismen en waarbij uiteindelijk vier zaadcellen met ieder de helft van het oorspronkelijke DNA-materiaal worden gevormd, is voor dit project irrelevant en zal dus ook niet behandeld worden. Wel gaan we in op het begrip apoptose, het principe van geprogrammeerde celsterfte. Hierbij worden levende cellen die overbodig of gemuteerd zijn, geforceerd om zichzelf te vernietigen.

2.4.1

Celcyclus

Zoals hiernaast schematisch is weergegeven, is de cyclus grofweg onder te verdelen in twee verschillende fases: ˆ De interfase waarin de cel groeit. Deze is weer onder te verdelen in drie fasen:

– G1: de eerste tussenfase, waarin de cel extra cytoplasma en eiwitten produceert; – S: de synthese-fase, waarin de chromosomen gerepliceerd worden. Hierdoor zullen na deling twee -op mutatie na- genetisch identieke cellen ontstaan; – G2: de cel bevat nu tweemaal zoveel genetisch materiaal en bereidt zich voor op de celdeling door de laatste benodigde stoffen te produceren.

Figuur 2.4: De celcyclus [4]

ˆ De mitotische fase of M-fase waarin de cel zich deelt. Deze bestaat uit:

– de mitose waarin de gekopieerde strengen DNA zich van elkaar scheiden, en – de cytokinese waarin de cel zich splitst in twee nieuwe cellen. Tussen de G1- en S-fase en tussen de G2- en M-fase bevinden zich twee zogenaamde checkpoints, respectievelijk C1 en C2. Bij C1 wordt bepaald of de betreffende cel de synthese-fase ingaat of niet. Als het DNA bijvoorbeeld niet meer intact is kan de cel vernietigd worden door apoptose. Als deling van de cel nog niet gewenst is, kan deze een tijd lang in G1 blijven. Ook kan de cel uit de celcyclus “stappen” en in een passieve status G0 komen waarin geen groei of deling meer plaatsvindt en waarin de cel zich gedurende lange tijd kan bevinden. Sommige cellen, zoals bijvoorbeeld zenuwcellen, verblijven voor de duur van het leven van het organisme in deze G0-fase, terwijl epitheelcellen in de darmen meer dan twee maal per dag delen. Het verschil in duur van celcycli van cellen komt vooral door de variatie van de tijdsduur die een cel doorbrengt in G1 en G0. Als de cel eenmaal voorbij checkpoint C1 is, zal de rest van de celcyclus zich relatief snel afspelen - bij zoogdieren binnen 12-24 uur[5]. Bij 6

Tumorvorming en kanker het punt C1 krijgt de cel dus het “startsignaal” om te gaan delen en een gehele celcyclus te volmaken. Dit zal van belang zijn voor het wiskundige model dat later wordt ge¨ıntroduceerd. Bij checkpoint C2 wordt gekeken of in de S-fase en de G2-fase al het DNA gerepliceerd is en of dit DNA nog intact is. Als dit zo is, is de cel klaar voor deling en zal hij over gaan naar de M-fase. Belangrijk is op te merken dat de cel groeit gedurende heel de interfase en tijdens de M-fase niet meer. Bovendien kan een cel in de M-fase ook niet migreren.

2.4.2

Apoptose

De cellen van meercellige organismen worden volgens complexe systemen aangezet tot delen of het stoppen daarmee. Echter, niet alleen het aantal cellen dat zich zal delen wordt gereguleerd, ook de celsterfte wordt gestuurd. Als cellen niet langer nodig zijn, kan een systeem van ‘geprogrammeerde zelfmoord’ geactiveerd worden. Dit heet apoptose 1 . De hoeveelheid cellen die door apoptose te gronde gaan is hoog. In een gezond volwassen mensenlichaam gaan elk uur miljarden cellen dood. Toch blijft in gezond volgroeid weefsel het aantal cellen min of meer constant en is het aantal cellen dat door apoptose sterft in balans met het aantal cellen dat gevormd wordt door mitose. Als een weefsel door externe factoren, bijvoorbeeld ziekte, opeens in omvang toeneemt, zal apoptose er uiteindelijk voor zorgen dat dit weer naar het oude niveau wordt teruggebracht. Andersom zal na verlies van weefsel de celdeling de overhand krijgen boven de celsterfte en De functie van de hoge frequentie waarmee cellen zich delen is het repareren van het lichaam en het hernieuwen van delen hiervan. We betalen hier echter ook een prijs voor.

2.5

Tumorvorming en kanker

Om het hele proces van celdeling en -sterfte in balans te houden, moet de individuele cel zich gedragen naar de behoeften van het gehele organisme. De cel moet sterven of delen, naargelang nieuw weefsel nodig is of oud weefsel juist afgebroken moet worden. Als beide niet het geval zijn, moet de cel in passieve stand zijn; hij moet zijn functie blijven uitvoeren en blijven op de locatie waar hij nodig is. Soms gehoorzaamt een cel niet aan hetgeen hem opgedragen is en functioneert hij niet in het belang van het gehele organisme. Het lichaam bestaat uit z´o veel cellen en is zo vindingrijk dat het geen significant kwaad ondervindt als dit bij ´e´en individuele cel gebeurt. Grote problemen en chaos in het lichaam kunnen echter ontstaan, als dit het gevolg is van een verandering in het DNA van een cel. Dit leidt tot ongecontroleerd delen, waarbij de dochtercellen dit genetisch mankement overerven en zich even ‘asociaal’ gaan gedragen. In het weefsel ontstaan dan tumoren: samengeklusterde cellen die zich snel delen. Als deze cellen op hun plek blijven, spreken we van een goedaardige tumor. Als ze zich echter uitbreiden naar omliggend weefsel en hier schade aanbrengen is de tumor kwaadaardig en lijdt het organisme aan kanker. Als deze cellen zich naar andere plaatsen in het lichaam verplaatsen en daar nieuwe tumoren vormen spreken we van een metastase of uitzaaiing. Een gezonde cel verandert pas in een kankercel als een reeks van mutaties zijn opgetreden in meerdere genen die van belang zijn bij het reguleren van de celdeling. De belangrijkste van deze genen zijn: 1

Het woord apoptose stamt van het Griekse ‘apo ptosis’, wat ’afvallen’ (als bij bladeren van een boom) betekent.

7

Inleiding in de Celbiologie ˆ Proto-oncogenen zijn betrokken bij het stimuleren van normale celdelingen. Ze zetten hun eigen cel of andere cellen aan tot groeien en uiteindelijk tot delen. Na mutatie kunnen deze genen veranderen in oncogenen en zetten ze aan tot tot overmatige deling of zelfs onbeperkte groei. ˆ Tumorsuppressorgenen zijn - zoals de naam al doet vermoeden - genen die er voor zorgen dat cellen niet buitensporig hard gaan delen en tumoren vormen. Mutatie kan er voor zorgen dat tumorvorming minder effectief wordt tegengegaan. ˆ Apoptose-genen regelen de geprogrammeerde celsterfte. Na mutatie kan dit zelfvernietigingsmechanisme haperen waardoor de kankercellen niet vernietigd worden. ˆ DNA-herstel regelende genen zorgen dat bovenstaande en andere mutaties ongedaan kunnen worden gemaakt. Als dit niet meer werkt kunnen steeds meer defecten in het DNA ontstaan.

Tumoren en kanker kunnen in bijna heel het lichaam ontstaan. Veel voorkomend zijn tumoren in de darmen, de longen of (bij mannen) in het prostaat en (bij vrouwen) in de borst.

2.5.1

Oorzaak

Kanker is een ziekte waarbij de oorzaak ligt in een opeenvolging van veranderingen - mutaties - in het DNA van een individuele cel. Bij het delen van cellen wordt DNA gekopieerd, maar soms gaat hier iets fout en wordt verkeerde informatie ‘overge¨erfd’. In veel gevallen heeft dit geen ernstige gevolgen- soms kan een mutatie zelfs een voordeel voor het organisme betekenen-, maar als dit net een aantal maal op de verkeerde plaats in het DNA gebeurt, kan een cel ongecontroleerd snel gaan delen en een tumor vormen. De grootste oorzaak van kanker ligt dus in het niet 100% accurate mechanisme van DNA-replicatie - overigens hetzelfde mechanisme dat mede bepaald heeft dat evolutie van eenvoudig naar complex organisme heeft plaatsgevonden en dat heeft gezorgd voor een enorme biodiversiteit. Het is dus een bijproduct van de manier waarop leven gevolueerd is en daarom zal tumorvorming waarschijnlijk nooit helemaal verdreven kunnen worden. Er zijn een hoop externe factoren die dit proces van gen-mutaties kunnen be¨ınvloeden. Een enkele vorm van kanker wordt deels veroorzaakt door een virus, maar voor de meeste geldt dit niet: kanker is geen infectieziekte. Statistisch epidemiologisch onderzoek [11] heeft aangetoond dat de omgeving een veel grotere rol speelt bij het ontstaan van kanker. Zo verschilt per land welke vorm van kanker vaker of minder vaak voorkomt. Dit hangt samen met de leefgewoonten en de natuurlijke omgeving van dit bepaalde gebied. Het lijden aan overgewicht bijvoorbeeld, is positief gecorreleerd met een hogere kans op verschillende kankers en de radioactieve straling in de zon verhoogt de kans op huidkanker. De grootste oorzaak van het ontstaan van kanker echter, is het roken van tabak. Geschat is dat het aantal kanker-gerelateerde sterfgevallen met 30% zou kunnen worden teruggedrongen als niemand meer zou roken[5]. Tabak zou een rol spelen bij de ontwikkeling van meerdere vormen van kanker, maar is bovenal verantwoordelijk voor bijna alle gevallen (85%) van kanker in de longen[12]. Volgens de meest recente cijfers [14] zorgde kanker in 2008 wereldwijd voor 7,9 miljoen sterfgevallen, 13,1% van het totaal. In datzelfde jaar zijn wereldwijd 12,7 miljoen nieuwe kankergevallen ontstaan, waarvan 6,6 miljoen mannelijk en 6,1 miljoen vrouwelijke slachtoffers[13]. Longkanker komt het meest voor en is met 1,6 miljoen goed voor 13% van de totale nieuwe kankergevallen. 8

Tumorvorming en kanker

2.5.2

Longkanker

Als er een kwaadaardige tumor ontstaat in weefsel van de longen, spreken we van longkanker. Dit is een verzamelnaam voor een aantal types die verschillen in locatie in de long en het type cel dat het betreft. De meeste hiervan zijn carcinomae, kwaadaardige tumoren die ontstaan in de epitheelcellen in de longen. Hierbinnen is weer een onderverdeling in celgrootte en het type epitheelcel, maar wij gaan ons in het onderzoek toeleggen op plaveiselcelcarcinoom, een kwaadaardige tumor die ontstaat vanuit plaveiselepitheel in de longen. Dit type longkanker treft zo’n 14.4 per 100.000 mensen per jaar [12] en vormt daarmee 30% van alle longkanker. In de jaren ’50 is vastgesteld [16] dat actief roken de belangrijkste oorzaak is van longkanker. Slechts in 15% van de gevallen van longkanker is het slachtoffer een niet-roker. In tabaksrook zitten zo’n 60 verschillende kankerverwekkende stoffen. Daarnaast onderdrukt nicotine het afweersysteem dat de ontwikkeling van tumoren tegengaat. Andere oorzaken van longkanker kunnen het passief meeroken en het inademen van asbest zijn. Plaveiselcelcarcinoma zijn in verhouding met andere vormen van longkanker vaker het gevolg van het roken van tabak, maar ook het humaan papillomavirus kan een rol spelen. Belangrijke symptomen van longkanker zijn het ophoesten van bloed, gewichtsverlies en kortademigheid. Ook ontstekingen in de longen kunnen optreden. De diagnose wordt gesteld op basis van radiografie of CT-scan en een biopsie. De behandeling bestaat uit een combinatie van chemotherapie, radiotherapie en chirurgie waarbij een longkwab of de gehele long wordt verwijderd. Longkanker is de meest voorkomende oorzaak van kankergerelateerd overlijden. Na het stellen van de diagnose overlijdt zo’n 60% van de pati¨enten binnen een jaar. De vijfjaars-overleving van alle soorten longkanker samen ligt op minder dan 10%[15].

2.5.3

De tumorcel

De delingseigenschappen van tumorcellen verschillen in grote mate van die van gezonde cellen. Ook zijn ze anders van vorm en de loopt de intercellulaire communicatie en de migratie van kankercellen niet zoals het hoort. Dit alles is essentieel voor de aard van de kanker. Hieronder volgt een beschrijving van de manier waarop tumorcellen afwijken van gezonde cellen. Morphologische verschillen tussen epitheel- en tumorcellen Gezonde plaveisel-epitheelcellen gaan allemaal min of meer hetzelfde groeiproces door. Hun grootte zal dan ook binnen de fysiologische grenzen van dit proces vallen en grote afwijkingen komen zelden voor. Kankercellen groeien zeer ongecontroleerd en dit kan men ook terugzien in de grootte en vorm van de cellen, die erg varieert. Biomechanica en biofysica van tumorcellen De stijfheid en vervormbaarheid - de mate waarin de cel vervorming toelaat als een kracht wordt uitgeoefend - zijn ook eigenschappen waarin kankercellen kunnen afwijken van gezonde cellen. Of dit in positieve of negatieve mate gebeurt hangt af van het type cel, zo blijkt uit onderzoek[18]. Een arts kan soms de tumor met de hand voelen doordat de stijfheid van het weefsel is veranderd. In onderzoek[19] is voor een longcarcinoom bij muizen vastgesteld dat de stof stikstofmonoxide (NO)- dat een rol speelt bij het vormen van metastasen- er voor zorgt dat deze cellen minder deformeren. Deze cellen zijn dus een stuk stijver dan gezonde cellen. 9

Inleiding in de Celbiologie De celcyclus van tumorcellen Zoals hierboven vermeld, verandert een gezonde cel pas in een kankercel na een serie mutaties, waarbij verschillende systemen die de celdeling regelen anders gaan functioneren dan eigenlijk zou moeten. Hierdoor wijkt de celcyclus van een kankercel op meerdere manieren af van die van een gezonde cel. Doordat tumorsuppressorgenen muteren zullen cellen bij checkpoint C1 niet meer tegengehouden worden als DNA defect is of celdeling niet gewenst is. Cellen zullen ook vaker tot delen overgaan door de mutatie van oncogenen. Ook worden de cellen niet meer door apoptose vernietigd door deze oncogenen en de gemuteerde apoptose-regulerende genen. Dit samen zorgt er voor dat de G1-fase van kankercellen korter is en ze minder snel vernietigd worden. Gezonde cellen sterven doorgaans net zo vaak als dat er nieuwe cellen bijkomen. Het aantal gezonde cellen zou onder normale omstandigheden dus gelijk blijven. Voor kankercellen geldt dit niet. Netto komen er meer cellen bij dan er sterven waardoor de kolonie zich uitbreidt, uiteindelijk ten koste van gezonde cellen. De verdubbelingstijd - de tijd waarin een kolonie zich heeft verdubbeld - van een kolonie tumorcellen ligt tussen de 100 en 200 dagen en gedurende deze tijd wordt het volume van deze tumor dus ongeveer 2 keer zo groot [22]. Doordat de tumor groeit zullen gezonde cellen door apoptose en verdrukking sterven en zal hun aantal afnemen. Intercellulaire communicatie tussen tumorcellen Gezonde cellen communiceren met elkaar via chemische en mechanische signalen. De wand van deze cellen is zeer permeabel, waardoor ionen en moleculen makkelijk van de ene cel in de andere kan overvloeien. Hierdoor kunnen ze elkaar aansporen tot bewegen als dit nodig is, tot groeien en delen als het weefsel beschadigd is en tot apoptose als er weefsel vernietigd moet worden of de cel defect is. De wanden van tumorcellen zijn relatief impermeabel[20] en bieden dus in veel mindere mate de mogelijkheid om deze chemische signalen te versturen en ontvangen. De intercellulaire communicatie die de groei stuurt functioneert nauwelijks en ook migratie van de tumorcellen wordt niet be¨ınvloed door de omliggende cellen. Zowel wat betreft de deling als wat betreft beweging trekken de kankercellen zich dus weinig aan van wat er om hen heen gebeurt.

2.5.4

Tumorgrootte

Voor dit onderzoek is het relevant om te weten wat de orde van grootte van een tumor is en een idee te hebben hoeveel tumorcellen dit zijn. Een tumor wordt doorgaans pas ontdekt als hij groter is dan een centimeter in doorsnede. Dit is het geval als er een miljard tumorcellen aanwezig zijn. Dit is bij een verdubbelingstijd van 150 dagen het geval na ongeveer log2 (1.000.000150 ) = 150 log2 (106 ) = 900 · 3.32 ≈ 2990 dagen ≈ 8.1 jaar. Cellen in de longen hebben door de ruime aanwezigheid van zuurstof een verdubbelingstijd van ongeveer 100 dagen. Doordat het weefsel heel dun is en al uitzaaiingen plaats kunnen vinden vanaf 50.000 cellen is longkanker dan ook een van de meeste gevaarlijke kankers. Zo’n uitzaaiing kan al plaatsvinden na 4,3 jaar.

10

Hoofdstuk 3

Mathematisch model Om het ontstaan van een tumor binnen een kolonie van lichaamscellen te kunnen simuleren, wordt een wiskundig model gemaakt. Om tot zo’n model te komen zullen eerst een aantal aannames omtrent de individuele cel en de kolonie als geheel moeten worden gedaan. Wiskundige formules beschrijven vervolgens hoe deze cellen onderling met elkaar communiceren en zo de migratie van cellen binnen een kolonie kunnen sturen en hoe cellen bewegen onder invloed van elkaar. Verder wordt gemodelleerd hoe cellen groeien, delen en sterven en hoe tumorcellen zich gedragen en uit kunnen groeien tot een kankergezwel. Hieronder komen de gedane modelaannamen aan bod en wordt de wiskundige uitwerking van het model nader toegelicht.

3.1

Epitheelweefsel in de longen

In dit project is gekeken naar een specifiek weefsel en specifieke cellen in het menselijk lichaam, namelijk de epitheelcellen in het oppervlakteweefsel van de longen. Deze situatie leidt tot bepaalde aannames op cel- en kolonieniveau die hieronder eerst zullen worden besproken. Leefomgeving Het celweefsel dat in het model wordt beschreven is het longepitheel; de dunne laag cellen die het oppervlak van de longen vormen dat grenst aan de lucht die ingeademd wordt. Deze laag is ´e´en cel dik. We zullen dit als een twee-dimensionaal oppervlak in R3 benaderen. De rand van het domein waarop de beweging van de cellen plaatsvindt is een cirkel. De omgeving waar de cellen zich in bevinden is dus ook te interpreteren als een laboratorium-experiment, waarbij de cellen naast elkaar op een voedingsbodem op een petrischaaltje liggen. Het middelpunt van deze leefomgeving is de oorsprong, (0, 0). De straal van het leefgebied van de kolonie is Rkol . Het domein waarbinnen de cellen zich kunnen bewegen en de kolonie zich uit kan breiden is dan: Ω = {(x, y) ∈ R2 : x2 + y 2 < Rkol } Het substraat waarop de cellen zich bevinden is in dit model homogeen, met andere woorden: het heeft op iedere locatie dezelfde samenstelling en dus ook dezelfde fysische eigenschappen.

Mathematisch model De cel De kolonie bestaat uit n cellen, waarbij het getal n varieert door de tijd als gevolg van deling en sterfte. De locatie van het middelpunt van cel i wordt op tijdstip t gegeven door ri (t) = (xi (t), yi (t)) ∈ Ω,

met i ∈ {1, 2, ..., n}.

De wand van de cellen wordt homogeen verondersteld: de cellen hebben allemaal dezelfde stijfheid en reageren op dezelfde manier op omliggende cellen. Morfologie van de cel Cellen binnen een weefsel kunnen onderling nog erg verschillen in grootte en vorm. Bovendien vervormt een cel bij migratie en deukt hij in door omliggende cellen. Er is dus niet ´e´en standaard vorm of grootte. Met een microscoop gemaakte beelden (zie figuur 3.1) van een cel in neutrale, niet-actieve toestand, tonen de projectie van een cel op het platte vlak. Dit laat een ronde vorm zien en de cel wordt dan ook als zodanig op de platte ondergrond verondersteld. In het model is de vorm van de cellen daarom genomen als zijnde een halve bol die met de platte kant op de ondergrond - het substraat- ligt. Het indeuken van de celwand door omliggende cellen of de rand van het domein vindt wel plaats, maar hier wordt bij de locatie van de cel niet gecorrigeerd voor het veranderde middelpunt.

Figuur 3.1: Kolonie van fibroblasten. Bij A ziet men de ovaal uitgerekte vorm van een migrerende cel. Bij B liggen de cellen in een dicht rooster en kan men de “ronde” vorm zien. [21]

De individuele cel i wordt dus weergegeven als een cirkel rond punt ri met straal Ri . In het model groeien de cellen. Aan de afmetingen is wel een grens: R0 ≤ Ri ≤ Rmax , waarbij R0 de straal is van een dochter cel gelijk nadat deling heeft plaatsgevonden, en Rmax de straal is van de moedercel vlak voor deling. Behoud van totaal volume voor en na deling geeft: 1 4 3 1 4 · πRmax = 2 · · πR03 2 3 2 3

3.2

⇒

Rmax =

√ 3

2R0 ≈ 1.26R0 .

Beweging van cellen

Zoals in de vorige sectie is behandeld, communiceren cellen met elkaar via trekkracht op het substraat en leidt dit tot het actief migreren van cellen. Als cellen door deze actieve migratie dusdanig naar elkaar toe zijn bewogen dat ze botsen, oefenen ze een tegengestelde kracht uit. De krachten die hierbij een rol spelen en de bewegingsvergelijkingen die dit beschrijven zullen hieronder kort toegelicht en uitgewerkt worden. Aanvankelijk wordt de situatie beschouwd waarin cellen nog niet tegen elkaar aan botsen. 12

Beweging van cellen

3.2.1

Bewegingsvergelijking

De beweging van een cel i wordt als volgt gegeven: dri (t) = Vi (t)b zi (t). dt

(3.1)

Hierbij is Vi (t) de absolute snelheid van cel i op tijdstip t en b zi (t) is de richting waarin de cel beweegt op dit tijdstip. De grootte van Vi en de richting van b zi (t) worden op ieder tijdstip bepaald door cel i zelf en de hieromheen gelegen cellen. Het wiskundig raamwerk dat deze beweging beschrijft wordt behandeld in de volgende paragrafen.

3.2.2

Actieve migratie

De beweging - migratie - van een cel bestaat uit 2 componenten, een actieve en een passieve. De actieve is de “kruipbeweging” die de cel maakt onder invloed van signalen van omliggende cellen. De passieve component is de beweging die wordt veroorzaakt door botsing met omliggende cellen of de rand van het weefsel. In deze Sectie zal de actieve migratie besproken worden. Rekenergie-dichtheid De trekkracht die een levende cel uitoefent, leidt tot kleine vervorming van de ondergrond. De rek zorgt ervoor dat potenti¨ele energie wordt opgeslagen in dit substraat. De energiedichtheid (energie per eenheid volume) is scalair en alle bijdrages geleverd door de verschillende cellen kunnen dan ook opgeteld worden. De energiedichtheid op de locatie van de cel bepaalt hoe snel hij zich verplaatst en de gradi¨ent van dit veld bepaalt in welke richting hij zich verplaatst. De door cel j ∈ {1, 2, ..., n} uitgeoefende trekkracht wordt in het model beschouwd alsof hij geleverd wordt door een puntbron in rj , het centrum van de cirkelvormige projectie van de cel op het substraat. De bijdrage die deze trekkracht levert aan de totale potenti¨ele energie op dit punt wordt gegeven door 1 Mj0 = σ, 2 0 waarbij Mj de potenti¨ele energie dichtheid is in cel-centrum rj ten gevolge van de door cel j uitgeoefende trekkracht Fj en σ en  respectievelijk de spanning en rek in het celcentrum zijn. Voor trekkracht Fj wordt evenredigheid met de oppervlakte van de projectie van de cel op het substraat aangenomen. Dit komt neer op: ( R Fb( R0j )2 als cel j levend, Fj = 0 als cel j dood. Hierbij zijn Fb en R0 respectievelijk de trekkracht en de radius van een cel gelijk na deling. Verder wordt aangenomen dat de elasticiteit lineair is. Voor kleine vervormingen van het substraat geeft de wet van Hooke dan: σ = Es (rj ), waarbij Es de elasticiteitsmodulus van het substraat is op locatie rj . Bovendien geldt nu voor de spanning: σ=

Fj Fj = , Aj πRj2 13

Mathematisch model met Aj de oppervlakte van de projectie van cel j. Bovenstaande vergelijkingen samenvoegen geeft nu voor een levende cel j: 1 Fb2 . (3.2) Mj0 = 2 2π Es R04 De elasticiteitsmodulus van het substraat is constant genomen wegens homogeniteit. Uit vergelijking (3.2) volgt dat Mj0 onafhankelijk is van de grootte van cel j

Uitdoving signaal over het substraat Het substraat wordt maximaal vervormd in het centrum van de cel. Verder van dit centrum af wordt deze vervorming steeds minder en dooft de dichtheid van de potenti¨ele energie en dus de sterkte van het het mechanisch signaal - uit. Een analytische oplossing van deze uitdoving in termen van Bessel-functies wordt hier benaderd door een exponentieel verband. De potenti¨ele energiedichtheid M (r) op plaats r is een som van alle bijdrages van omliggende cellen op dat punt. De bijdrage van cel j wordt dan benaderd door Mj (r) = Mj0 exp{−λj

||r − rj || }, Rj

s met λj = E Ec , de verhouding tussen de elasticiteitsmoduli van respectievelijk het substraat en de cel. Homogeniteit geeft λj = λ. Een plot hiervan is te zien in figuur 3.2:

Figuur 3.2: Uitdemping trekenergie-dichtheid In het model is voor actieve migratie van cel i alleen van belang wat er gebeurt op locatie ri . Dit geeft voor cel i in een kolonie van n cellen: M (ri ) =

n X j=1

Mj (ri ) =

n X

Mj0 exp{−λ

j=1

||ri − rj || }. R

In figuur 3.2 is duidelijk te zien dat het mechanisch signaal zeer snel uitdooft en dus Mi (ri ) >> Mj (ri ), voor alle j 6= i. Voor twee aan elkaar grenzende cellen j en i is deze verhouding als volgt: ||R +R ||

Mj0 exp{−λ jRj i } 2R Mj (ri ) −λ R 0 − λ max = e 1.26 = ≤ e ||0|| 0 Mi (ri ) Mi exp{−λ Ri }

(≈ 1.28 · 10−7 met λ = 10). (3.3)

Dit heeft tot gevolg dat de potentiaal op punt ri nagenoeg alleen wordt bepaald door de trekkracht van cel i zelf, dus M (ri ) ≈ Mi (ri ) = Mj0 . 14

Beweging van cellen Hieruit volgt dat de snelheid waarmee de cel zich voortbeweegt - indien niet gehinderd door omliggende cellen- nagenoeg constant is. Detectie drempel Signalen van cellen die ver van elkaar af liggen kunnen zo zwak zijn dat de cel deze niet meer opvangt. Er geldt dan ook een drempel voor dit signaal. De afstand waarover epitheel-cellen elkaars aanwezigheid nog kunnen voelen hangt af van de stijfheid van het substraat. Uit studie van Reinhardt-King [9] blijkt dat deze afstand voor een vergelijkbaar substraat rond de dˆ = 30µm ligt. Dit komt voor het longepitheel overeen met een minimaal detecteerbare energiedichtheid van dˆ ε = Mi0 exp{−λ } ≈ 2.7 · 10−31 R Dit geeft voor het signaal dat cel i voelt ten gevolge van de trekkracht van cel j de volgende vergelijking: ||ri − rj || }1||ri −rj ||≤30µm R de indicatorfunctie zodat  0 als ||ri − rj || > 30µm 1||ri −rj ||≤30µm = 1 als ||ri − rj || ≤ 30µm Mj (ri ) = Mj0 exp{−λ

Hier is

3.2.3

1||ri −rj ||≤30µm

Passieve migratie

Cellen trekken elkaar niet alleen aan wanneer ze bij elkaar in de buurt liggen maar stoten elkaar ook af als ze met elkaar botsen. Het indeuken van de celwand zorgt voor een afstotende kracht die het naar elkaar toe migreren van de cellen tegenwerkt. Deze kracht werkt dus negatief op de bijdrage aan potenti¨ele energiedichtheid van de twee cellen op elkaar.

Figuur 3.3: Botsende cellen i en j Voor cel i en j, i 6= j, die elkaar een afstand h = grootte van deze bijdrage gegeven door [23]: M

ij

4Ec = √ 15 2π

1 2

max(0, Ri + Rj − |ri − rj |) indeuken is de



h Ri

5

2

(3.4) 15

Mathematisch model Deze potenti¨ele energie dichtheid moet nu afgetrokken worden van de positieve bijdrage ten gevolge van het aantrekken van de twee cellen. Dit geeft voor de totale rek-energiedichtheid die werkt op cel i: M (ri ) = Mi0 +

n X

|Mj (ri ) − M ij |.

(3.5)

j=1,j6=i

Figuur 3.4: M ij als gevolg van bostende cellen i en j.

Weefselwand De rand van het weefsel wordt beschouwd alsof hij bestaat uit hetzelfde materiaal als de celwand zelf. Het tegen de wand aan “botsen” geeft dezelfde reactie als tegen een naburige cel aanbotsen. Omdat is aangenomen dat de wand van het weefsel een perfecte cirkel met straal Rkol is, wordt de indeuking van een cel door contact met de wand gegeven door: hw =

1 max(0, ||ri || + Ri − Rkol ). 2

Analoog aan het botsen met een naburige cel veroorzaakt dit een bijdrage aan de rekenergiedichtheid:  5 4Ec hw 2 i . Mwand = √ 15 2π Ri

3.2.4

Mitotische cellen

Cellen die zich in de mitotische fase bevinden en zich voorbereiden op de deling, voeren geen actieve migratie meer uit. Ze bewegen nog wel onder invloed van botsende buurcellen en botsingen met de celwand, maar niet als gevolg van de rek-energiedichtheid afkomstig van omliggende cellen. Ondanks dat cel j het substraat wel degelijk vervormt op plaats ri , is vanuit het oogpunt van cel i de rek-energiedichtheid op deze locatie dus gelijk aan nul. De totale potenti¨ele energie-dichtheid op locatie ri wordt nu hergedefini¨eerd als de effectieve energie-dichtheid, met andere woorden: de potenti¨ele energie-dichtheid die ook daadwerkelijk beweging van cel i tot gevolg heeft. Dit geeft voor vergelijking (3.5): M (ri ) = Mi0 +

n X j=1,j6=i

16

|ψi Mj (ri ) − M ij |,

(3.6)

Beweging van cellen waarbij  ψi =

3.2.5

0 1

als cel i mitotisch als cel i niet mitotisch.

Bewegingsrichting

De bewegingsrichting van de cel wordt over een tijdsinterval ∆t gegeven door een vector. Deze vector is een lineaire combinatie van alle eenheidsvectoren die wijzen in de richting van de omliggende cellen. De eenheidsvector van cel i naar cel j (met i 6= j) wordt gegeven door: vij =

rj − ri . ||rj − ri ||

In het geval van cellen die in contact zijn met de wand van het cirkelvormig veronderstelde weefsel is deze vector normaal aan de wand: vi,w = −

ri . ||ri ||

De netto bewegingsrichting van de cel wordt nu bepaald door een gewogen gemiddelde van al deze eenheidsvectoren. Het gewicht dat wordt toegekend aan ieder van deze vectoren is gelijk aan de totale potenti¨ele energie die cel i ondervindt van cel j. Dit geeft voor de richting van beweging van cel i de volgende formule: i vi,w + zi = Mwand

X


vij ψi Mj (ri ) − M ij .

i6=j

De bijdrage is negatief als ψi Mj (ri )−M ij < 0 en de bewegingsrichting is dan ook tegengesteld aan vij . De uiteindelijke richting van bewegen wordt nu gegeven door de genormaliseerde vector: zi zbi = . (3.7) ||zi ||

3.2.6

Beweging numeriek ge¨ıntegreerd

Bovenstaande formules voor de rek-energie (3.5) en de bewegingsrichting (3.7) leiden nu tot een discretisering van bewegingsvergelijking (3.1). De snelheid waarmee cel i zich beweegt op tijdstip t wordt gegeven door [1]: Vi = αi M (ri ), 3

waarbij parameter αi [ m N s ] geschreven kan worden als αi =

βi Ri3 Fi , cf Fb2

met mobiliteitsparameter βi [s−1 ] en wrijvingsco¨effici¨ent cf . Er wordt aangenomen dat Fi = Fb. De beweging van cel i ∈ {1, ..., n} over tijdsinterval ∆t wordt nu gegeven door de volgende numerieke integratie: R3 (3.8) ri (t + ∆t) = ri (t) + ∆tβi i M (ri )zbi , cf Fb 17

Mathematisch model

3.3

Cel cyclus

Bovenstaand bewegingsmodel is grotendeels gebaseerd op het semi-stochastische cel-gecentreerde model voor migratie zoals dat beschreven is in [1]. In dit oorspronkelijke model zijn de celdeling en -sterfte als volledig stochastisch verondersteld; op ieder tijdstip heeft elke cel in de kolonie een bepaalde kans om te sterven en om te delen. De cellen maken hierbij niet de cyclus van synthese en mitose door. In dit project is gekeken hoe dit op een meer realistische manier kan worden gemodelleerd en op een effici¨ente manier numeriek kan worden ge¨ımplementeerd. Hierbij worden de groei van de cellen, de verschillende fases en de celdeling meegenomen en verder worden deze gekoppeld aan de druk die de cellen ondervinden van omliggende cellen. In het model is zo goed mogelijk geprobeerd de verschillende fases van de cyclus mee te nemen en te beschrijven hoe de lichaamscel zich gedraagt in ieder van die toestanden. De cel kent in het model vier verschillende toestanden φi , waarvan drie in de interfase en een in de mitotische fase. 1. G1 . De eerste groeifase waarin een nieuwe dochtercel zich onmiddellijk na deling bevindt. In deze fase groeit de straal van de cel van R0 naar RC . Deze groei gaat stochastisch met lineaire trend. De cel migreert actief en kan sterven. De gemiddelde tijd waarin de cel zich in deze fase bevindt is TG1 . 2. S. De rustfase bij checkpoint 1. De cel wacht hier totdat hij kan beginnen aan de laatste groeifase. De cel neemt niet toe in omvang. Actieve migratie vindt wel plaats. De kans dat de cel overgaat in de volgende groeifase wordt be¨ınvloed door de druk die de cel ondervindt van omliggende cellen. De gemiddelde tijdsduur die een cel hierover doet is TS . 3. G2 . De tweede en laatste groeifase voordat de cel de mitotische fase ingaat. De straal groeit nu van RC naar RM en deze groei gaat stochastisch met dezelfde lineaire trend als in G1 . Het tijdsbestek waarover dit gebeurd is TG2 . 4. M. Dit is de mitotische fase waarin het DNA gescheiden wordt. De omvang van de cel neemt niet meer toe en aan het einde van deze fase deelt de cel zich tot 2 dochter cellen die zich beide weer in G1 bevinden en straal R0 hebben. Actieve migratie vindt niet plaats, de cel kan nog wel door omliggende cellen van haar plaats ‘gebotst’ worden. De cel bevindt zich gedurende een tijdsduur TM in deze mitotische fase alvorens zich te delen. Voor de verwachtte duur van de cyclus geldt nu TC = TG1 + TS + TG2 + TM . De cellen kunnen op ieder tijdstip sterven met een bepaalde kans die afhankelijk is van de druk van de omliggende cellen. De concentratie van cellen be¨ınvloedt dus zowel de delingssnelheid als het sterfteratio. Dit gaat een grote rol spelen als later in het project tumorcellen worden toegevoegd aan het model.

3.3.1

Groei

Om de groei te modelleren van de cellen, moet voor elke cel bepaald worden in welke fase van de cyclus hij zich bevindt. Voor φi = G1 en φi = G2 vindt toename in volume plaats, in S en M blijft het volume gelijk. De groeiparameter γi wordt gedefinieerd als volgt:  1 als φi = G1 ∨ φi = G2 γi = 0 als φi = S ∨ φi = M 18

Cel cyclus Vervolgens wordt hier aangenomen dat de straal van een groeiende cel met lineaire trend µ > 0 toeneemt. Dit geeft het deterministische deel van de differentiaalvergelijking: dRi = γi µ. dt Verondersteld wordt echter dat de groei van de cel niet deterministisch is en dat er dus ook een stochastische term in deze vergelijking zit. Om deze onzekerheid te modelleren, introduceren we eerst het principe van een Wiener Proces.

Wiener Proces[24] Een continu stochastisch proces W is een Wiener proces, ofwel random walk, als de volgende eigenschappen gelden voor W = W (t): ˆ W (0) = 0; ˆ Laat t0 < t1 < t2 , dan zijn gebeurtenissen W (t2 ) − W (t1 ) en W (t1 ) − W (t0 ) onafhankelijk; ˆ W (t) is stochastisch continu, m.a.w:

lim P (|W (t) − W (t0 )| > ) = 0,

t→t0

∀ > 0,

met P de bijbehorende kansmaat; ˆ Voor t 6= t0 , geldt W (t)−W (t0 ) ∼ N (0, t−t0 ) , m.a.w W (t)−W (t0 ) is normaal verdeeld met gemiddelde nul en variantie t − t0 .

We laten nu groei van Ri verstoord worden met variantie σ 2 t vanaf t0 = 0. Dan wordt de stochastische groei van de straal gemodelleerd door: dRi (t) = γi µdt + σdW (t).

(3.9)

Dit geeft voor de numerieke integratie met de voorwaartse methode van Euler: ∆Ri (t) = Ri (t) − Ri (t − ∆t) = γi µ∆t + σ (W (t) − W (t − ∆t)) . Hieruit volgt:   σ2 ∆Ri ∼ N γi µ, √ ∆t. ∆t De numerieke integratie van vergelijking (3.9) wordt dan gegeven door: σ2 ∆Ri = γi µ∆t + N (0, √ )∆t ∆t

(3.10)

met σ hier de standaardafwijking van de groei van een cel over ´e´en seconde.

19

Mathematisch model De verwachte groeisnelheid µ wordt in vergelijking (3.10) nu dusdanig gekozen dat de verwachtte totale tijd die het duurt dat cel i aan het groeien, is gelijk is aan TG1 + TG2 = TG . Dit geeft µ=

3.3.2

Rmax − R0 . TG1 + TG2

Checkpoint

Als de cel in fase G1 zit, groeit hij tot hij een bepaalde straal, RC , heeft bereikt. Wanneer dit het geval is, raakt hij in pauze-stand S waarin de cel een tijdje kan blijven alvorens verder te gaan met delen. Voor een cel in deze toestand geldt γi = 0, en dus: dRi (t) = σdW (t). Het verblijf van de cel in fase S wordt ook stochastisch verondersteld, en de kans waarmee een cel binnen tijdsinterval ∆t over kan springen naar G2 is Pg (∆t). Deze kans is dusdanig dat voor de verwachtte verblijfsduur TS,i van cel i in fase S geldt: E(TS,i ) = TS .

(3.11)

De kans dat cel i in fase φi = S binnen tijdsinterval [0, T ] overspringt naar fase G2 wordt in het discrete geval, T = n∆t gegeven door een geometrische verdeling: P (TS,i = n∆t) = (1 − Pg (∆t))n−1 Pg (∆t) Voor het verwachtte aantal tijdstappen n geldt dan: E(n) =

1 . Pg (∆t)

Dit geeft voor de verwachtte tijdsduur dat cel i in fase S blijft:

E (TS,i ) =

∞ X

Tn P (TS,i = n∆t) =

n=1

∞ X

n∆t(1 − Pg (∆t))n−1 Pg (∆t) = ∆tE(n) =

n=1

∆t . Pg (∆t)

Nu geeft conditie (3.11): Pg (∆t) = P (φi (t + ∆t) = G2 | φi (t) = S) =

∆t TS

Voor bepaalde keuzes van tijdsduren TG1 , TS en TG2 krijgen we de volgende plot:

Figuur 3.5: Groei van de straal van een gezonde cel en een kankercel. 20

(3.12)

Cel cyclus In figuur 3.5 is te zien hoe de straal van een gezonde cel gedurende een tijdsinterval van ongeveer 400 seconden groeit van R0 = 3 µm tot RC = 4 µm, waarna de straal zo’n 300 seconden gelijk blijft in grootte en tot slot in weer ongeveer 400 seconden groeit totdat zijn maximale grootte is aangenomen en de cel deelt tot 2 cellen, ieder met straal R0 = 3 µm. Bij de kankercel groeit de cel sneller en slaat hij bovendien de pauzestand over. Hij groeit dan ook in zo’n 300 seconden tot maximale straal en deelt dan tot 2 cellen met beginstraal R0 .

3.3.3

Mitose

Als cel i in fase G2 gegroeid is tot de maximale straal is aangenomen, Ri = Rmax , bevindt hij zich in de mitotische fase M. De cel gaat zich nu voorbereiden op de deling die aan het einde van deze fase plaatsvindt. De tijdsduur hiervan, TM , is een stuk kleiner dan die van de andere fases en wordt als deterministisch verondersteld. Gedurende tijdsinterval TM blijft de straal constant (γi = 0) op de stochastische term in vergelijking (3.10) na. Bovendien voert de cel geen actieve migratie uit; hij beweegt alleen onder invloed van botsende buurcellen. Voor de groei- en bewegingsparameters geeft dit voor mitotische cel i: γi = 0

en ψi = 0

Stel nu dat op tijdstip τ cel i zich gedurende een tijdsinterval TM in de mitotische fase heeft bevonden en klaar is voor deling. Cel i deelt zich tot dochtercellen ia en ib , waarvoor nu geldt dat de straal op tijdstip t = τ + dt gegeven wordt door Ri,a (τ + dt) = Ri,b (τ + dt) = R0 . De locatie waarop de moedercel zich zou bevinden in geval van geen deling, e ri (τ + dt), wordt het middelpunt tussen de locatie van de 2 dochtercellen, respectievelijk ri,a (τ + dt) en ri,b (τ + dt). De richting ubi waarin deling plaatsvindt, is uniform verdeeld. Dit geeft voor de locaties van de dochtercellen: ri,a (τ + dt) = e ri (τ + dt) + δR0 ubi = ri (τ ) + dtβi

R3 + δR0 ubi (τ + dt) cf Fi

ri,b (τ + dt) = e ri (τ + dt) − δR0 ubi = ri (τ ) + dtβi

R3 − δR0 ubi (τ + dt) cf Fi

met ubi (t) = (cos(2πξi ), sin(2πξi )), ξi ∼ U (0, 1), waarbij U (0, 1) de uniforme verdeling is op interval [0, 1].

Figuur 3.6: Cel i splitst zich op in cel i, a en cel i, b.

21

Mathematisch model

3.3.4

Apoptose

In het model z´ onder drukafhankelijkheid wordt de apoptose als volledig stochastisch verondersteld. De kans om gedurende tijdsinterval ∆t te sterven wordt dan aangegeven als Papt (∆t). Deze sterftekans wordt nu dusdanig gekozen, dat de kans dat een cel sterft gedurende tijdsinterval [0, T ], in de discretisatie onafhankelijk is van de grootte van tijdstap ∆t met T = n∆t. Voor de kans dat een cel tijdsinterval [0, T ] overleeft moet dan dus gelden T T · ∆t) = (1 − Papt (∆t)) ∆t = (1 − Papt (∆t))n . (3.13) ∆t Voor T = TC , de gemiddelde tijdsduur tussen twee delingen van “dezelfde” moedercel, geeft het omschrijven van vergelijking (3.13):

1 − Papt (n∆t) = 1 − Papt (

TC

(1 − Papt (∆t)) ∆t = 1 − Papt (TC ).

(3.14)

Hierbij is 1 − Papt (TC ) de kans dat een cel de gehele celcyclus doormaakt en zich deelt tot twee dochtercellen. We willen nu Papt (∆t) bepalen afhankelijk van de gemiddelde duur van een celcyclus (TC ) en de verdubbelingstijd van de kolonie (τ2 ). Stel nu dat op tijdstip t0 = 0 voor het aantal cellen n in de kolonie geldt: n(t0 ) = n0 , met n0 > 0. Dan geldt voor het verwachtte aantal cellen op tijdstip TC : E(n(TC ) | n(0) = n0 ) = 2(1 − Papt (TC ))n0 . (3.15) Voor Papt (TC ) <

1 2

wordt nu de verdubbelingstijd τ2 van de kolonie gedefinieerd: τ2 = {t ≥ TC : E(n(t0 + t)) = 2n(t0 ), ∀t0 ≥ 0}.

(3.16)

Samenvoegen van (3.15) en (3.16) geeft nu: τ2

E(n(τ2 ) | n(0) = n0 ) = (2(1 − Papt (TC ))) TC n0 = 2n0 .

(3.17)

Dit oplossen geeft: TC

1 − Papt (TC ) = 2 τ2

−1

en invullen van (3.14) geeft de gewenste sterftekans: Papt (∆t) = 1 − 2

( τ1 − T1 )∆t 2

C

.

(3.18)

In figuur 3.7 is de afhankelijkheid van de sterftekans Papt (TC ) en de verdubbelingstijd (blauwe lijn) en halveringstijd (rode lijn) van een celkolonie te zien. Voor een kolonie die in verwachting een constant aantal cellen houdt, dus Papt (TC ) = 0.5, geldt voor de verdubbelingstijd: τ2 → ∞. In de volgende paragraaf wordt de druk die de omliggende cellen uitoefenen ge¨ıntroduceerd. Het toenemen van deze druk zal een grotere kans op sterfte tot gevolg hebben. Dit zorgt voor een verandering in de exponent van vergelijking (3.18). 22

Figuur 3.7: Aantal celcycli voordat een cel verdubbelt of halveert

Invloeden van externe druk

3.4

Invloeden van externe druk

Om de vorming van een tumor en de dynamiek tussen gezonde cellen en tumorcellen goed te kunnen simuleren, moet eerst nog een tussenstap gemaakt worden: de aanwezigheid van omliggende cellen moet de duur van de celcyclus en de kans op apoptose be¨ınvloeden.

3.4.1

Celdruk

Om de druk die omliggende cellen uitoefenen op cel i te modelleren beschouwen we Paragraaf 3.2.3. Hier wordt de druk die cel j uitoefent op cel i gegeven door M ij , zie vergelijking (3.4). We defini¨eren nu de totale celdruk in cel i als Λi =

n X

i M ij + Mwand

(3.19)

j=1,j6=i

met als eenheid [ mN2 ] = [Pa].

3.4.2

Celcyclus

Zoals besproken in Sectie 2.4.1 wordt op de locatie van checkpoint C1 besloten of de cel verder gaat met groeien en zich voor gaat bereiden op mitose. Cellen in een kolonie vangen signalen op van omliggende cellen en reageren op de toestand in het weefsel. Als er een tekort is aan cellen, m.a.w er is ruimte om te delen, zullen cellen eerder geneigd zijn over te gaan naar de tweede groeifase, terwijl ze juist langer in passieve toestand S zullen blijven als de druk van omliggende cellen hoog is. In het hier besproken model komt dit overeen met een verlaging van kans Pg bij een toename van celdruk Λ. De statistische kwantificering van deze invloed is niet evident en daarom wordt in dit model een meer intu¨ıtieve benadering gebruikt. Gezocht wordt naar een dalende functie f (Λi ), zodat de kans dat cel i naar G2 overgaat gelijk is aan Pgi (∆t, Λi ) =

∆t f (Λi ). TS

Een functie die voldoet aan de wensen is f (Λi ) =

c1 . 1 + c2 expc3 Λi

Gekozen wordt nu om voor Pgi (0, Λi ) = 1.1 T∆tS te nemen, waarbij T∆tS als referentietoestand wordt gezien. De reden hiervoor is dat de deze referentietoestand moet gelden voor de gemiddelde cel en dat deze niet “drukloos” is. Dit geeft samen c1 = 11, c2 = 9 en c3 = 1. Pgi (∆t, Λi ) =

11∆t . TS (1 + 9 expΛi )

(3.20)

In figuur 3.8 is het verloop van de kans op groei afhankelijk van stapgrootte ∆t en de externe druk Λi te zien.

3.4.3

Drukafhankelijke sterfte

Voor de sterfte van een cel geldt nu de omgekeerde relatie: de kans op sterfte gedurende interval ∆t, Papt (∆t), neemt toe als de celdruk toeneemt. We zoeken nu een stijgende functie g(Λi ) zodat ( 1 − 1 )∆t · g(Λi ) i Papt (∆t, Λi ) = 1 − 2 τ2 TC . 23

Mathematisch model

Figuur 3.8: De kans om van fase S naar G2 over te gaan afhankelijk van de druk Deze functie wordt na proberen als volgt gekozen g(Λi ) =

9 1+10 exp−Λi

=

9 11 f (−Λi ).

Hierbij

i (∆t, 0) Papt

wordt de toestand gelijkgesteld aan 90% van de referentietoestand uit 3.18 en dit geeft voor de functie voor apoptose: i Papt (∆t, Λi ) = 1 − 2

( τ1 − T1 )∆t · g(Λi ) 2

C

.

Figuur 3.9: Drukafhankelijke sterftekans Papt (∆t) In figuur 3.10 is te zien hoe deze sterftekans afhangt van celdruk Λi .

3.5

Tumorvorming

Tumorcellen zijn anders van structuur, hebben een afwijkende celcyclus en bewegen op een andere manier dan gezonde cellen. Om de ontwikkeling van ´e´en tumorcel tot een tumor van enkele honderden cellen op een juiste manier te kunnen simuleren moeten deze eigenschappen worden ingebouwd in het model. Hieronder zal worden toegelicht hoe dit gebeurt. 24

Tumorvorming

3.5.1

De tumorcel

Eerst wordt nu in notatie onderscheid gemaakt tussen de gezonde cel en de tumorcel. Stel dat op tijdstip t er zich n(t) gezonde en n∗ (t) tumorcellen in de kolonie vinden. Dit maakt een totaal van N (t) = n(t) + n∗ (t) cellen. Iedere gezonde cel i ∈ {1...n(t)} bevindt zich op locatie ri (t) en iedere tumorcel k ∈ {1...n∗ (t)} bevindt zich op locatie r∗k (t). Ondanks dat kankercellen in werkelijkheid meer vari¨eren in vorm en grootte dan gezonde cellen, nemen we in dit model aan dat de kankercellen die op het substraat liggen ook de halve bolvorm hebben, waarvan de straal valt binnen: ∗ R0∗ ≤ Rk∗ ≤ Rmax .

3.5.2

Intercellulaire communicatie en migratie van tumorcellen

Tumorcellen worden minder be¨ınvloed door omliggende cellen en actieve migratie vormt een veel minder grote rol in het bewegingspatroon van cellen. In dit model wordt dan ook aangenomen dat de trek-energiedichtheid die wordt veroorzaakt door cellen die een trekkracht uitoefenen op het substraat niet wordt gedetecteerd door tumorcellen. Dit geeft vanuit het oogpunt van de tumorcel bekeken1 : Mj (r∗k ) = 0,

∀j ∈ {1, ..., N (t)}, k ∈ {1, ..., n∗ (t)}.

De beweging van de tumorcel wordt in dit model dus alleen bepaald door een botsing met buurcellen of de weefselwand. Hieruit volgt de beweging voor tumorcel k: N (t) j vj,w + zk = Mwand M (r∗k )b

X

−vk,j M kj .

(3.21)

j=1,j6=k

Met gehele bewegingsvergelijking dr∗k R3 = βi i M (r∗k )zbk dt cf Fk en met numerieke integratie r∗k (t + ∆t) = r∗k (t) + ∆tβk

Rk3 M (rk )zbk . cf Fk

(3.22)

met dezelfde mobiliteitsparameter βi [s−1 ] en frictieco¨effici¨ent cf als in vergelijking (3.8) voor de gezonde lichaamscel.

3.5.3

Celcyclus van een tumorcel

Zoals vermeld in Paragraaf 2.5.3 is de celcyclus van een tumorcel korter en de kans op sterfte kleiner. De groei gaat ongecontrolleerd en de invloed van omliggende cellen op de celgroei en de apoptose is dan ook veel kleiner. Het grootste verschil in de lengte van de celcyclus zit hem in de groeisnelheid ∗ ∗ TG1 + TG2 < TG2 + TG2 , 1 Nota bene: cel j heeft wel degelijk een bijdrage aan de energiedichtheid Mj (r∗k ) op locatie r∗k , maar omdat deze door tumor cel k niet wordt gedetecteerd wordt deze in het model beschouwd als afwezig. Ofwel, we stellen (Mi (r∗k ) = 0

25

Mathematisch model en het feit dat kwaadaardige tumorcellen volgens dit model geen ’pauzestand’ kennen: TS∗ = 0. Ook is de verdubbelingstijd van een kolonietumorcellen vele malen lager dan die van een kolonie gezonde cellen: τ2∗ < τ2 . Dit samenvoegen geeft voor de drukafhankelijke sterftekans:

Figuur 3.10: Drukafhankelijke sterftekans Papt (∆t) Als ook nog de verminderde gevoeligheid voor celdruk wordt meegenomen komt deze kans ook nog lager uit: ( 1∗ − 1∗ )dt · g ∗ (Λk ) ∗ Papt (dt, Λk ) = 1 − 2 τ2 TC (3.23) met g ∗ (Λk ) =

26

1.9 . 1 + e−3Λk

(3.24)

Hoofdstuk 4

Numerieke implementatie In dit onderzoek is de nadruk gelegd op het effici¨ent en biologisch realistisch implementeren van het ontstaan van een tumor. Dit is gedaan in het softwarepakket Matlab 2012b. De code met gedetailleerd commentaar is te vinden in bijlage A. Om de implementatie zo tijds-effici¨ent mogelijk te maken zijn alle berekeningen omgeschreven naar matrixnotatie. Matlab voert operaties op matrices vele male sneller uit dan wanneer dezelfde berekeningen als for-loops zijn ge¨ımplementeerd. In de uiteindelijke implementatie van het model is alleen de tijdsintegratie nog een for-loop.

4.1 4.1.1

Gebruikte Numerieke methode Beweging

Het uiteindelijke model voor de beweging van de cellen is een eerste orde differentiaalvergelijking in de tijd, zie vergelijking (3.1): r0i = Vib zi = αi M (ri )b zi . Vergeleken met het basismodel van Vermolen [1], is de stochastische term voor de snelheid hieruit weg genomen. In de eerste implementatie van het model voor tumorgroei
is gebruik gemaakt van een Euler voorwaartse methode in de tijd met een fout van O ∆t . Dit geeft de discretisatie: rn+1 = rni + ∆tαi M (rni )b zni . i

(4.1)

Voor de grote tijdstappen, bijvoorbeeld dt = 0.2, blijkt instabiliteit op te treden en de oplossingen te ’ontploffen’: de cellen leggen door de aantrekkende kracht en de grote tijstappen dusdanig grote afstanden af dat ze elkaar steeds meer indeuken en dus elkaar ook steeds verder wegstoten. De keuze ∆t = 0.1 geeft hier echter geen problemen en dus is verdere analyse van de stabiliteit niet nodig. Ook het variabel kiezen van de tijdstap, zoals gedaan wordt in vergelijking (3.1), blijkt na een aantal simulaties niet nodig te zijn. Om de implementatie tijdseffici¨ent te houden is gekozen om de voorwaartse methode niet te vervangen door een methode met hogere orde nauwkeurigheid of betere stabiliteit. Een impliciete methode zou eisen dat op ieder tijdstap een niet-lineair algebra¨ısch stelsel iteratief opgelost moet worden.

Numerieke implementatie

4.1.2

Deling

Als deling plaatsvindt, splitst cel i zich in twee nieuwe cellen: ia en ib . In Sectie 3.3.3 is dit al verder uitgewerkt. De discretisatie van deze vergelijking is een variant op vergelijking (4.1):

bi, rn+1 = rni + ∆tαi M (rni )b zni + δR0 u i,a bi. rn+1 = rni + ∆tαi M (rni )b zni − δR0 u i,b b is de eenheimet δ de factor waarmee beide cellen nog ingedeukt zijn op tijdstip t + ∆t. u dsvector in een random richting, uniform verdeeld over 360◦ .

4.1.3

Celgroei

De integratie van de groei van de cel wordt ook met
een Voorwaartse methode van Euler gedaan, waardoor de nauwkeurigheid hier ook O ∆t is. Zoals in Paragraaf 3.3.2 toegelicht bevat de celgroei een stochastische term. De numerieke integratie van deze groei wordt voor een niet-delende cel i gegeven door Rin+1 = Rin + γi µ∆t + σ(W n+1 − W n ),

(4.2)

met σ(W n+1 − W n ) ∼ N (0, σ 2 ∆t). Hiervoor wordt in de implementatie een randomizer gebruikt. Mitotische cel i geeft na deling twee dochtercellen ia en ib waarvoor geldt: n+1 n+1 = R0 . Ri,a = Ri,b

4.2

Implementatie

Hieronder volgt een korte toelichting tot de Matlab-implementatie van het model. De code is gedetailleerd becommentarieerd en te vinden in de Appendix. Om de berekeningen tijdseffici¨ent te houden zijn alle operaties in matrix-vector vorm gedaan.

4.2.1

Vectorisatie

Om hieronder de notatie wat te verduidelijken worden matrices in dikgedrukte hoofdletter (X) en vectoren met een dik gedrukte kleine letter (ri = (xi , yi )) of een hoofdletter met een streep eronder (Ri ) geschreven. Op tijdstip t wordt de locatie van het centrum van de cellen uit een kolonie van grootte N gegeven door locatiematrix: 

   x1 (t) y1 (t) r1 (t)    ..  .. .. X(t) =   =  . , . . xN (t) yN (t) rN (t) waarbij ri (t) = (xi , yi )(t) ∈ Ω de cartesische co¨ordinaten zijn van het centrum van cel i op tijdstip t. 28

Implementatie

4.2.2

Celcyclus

De manier waarop de cellen overgaan in de verschillende cycli is deels stochastisch en deels deterministisch. De fase waarin cellen {1, ...N } zich bevinden wordt gegeven door:  1 als cel i in fase G1    2 als cel i in fase S φk = {φ}i met φi = . 3 als cel i in fase G2    4 als cel i in fase M Stel nu dat op tijdstip tk voor cel i geldt: ˆ Rik < RC en φki = 1. Dan gaat cel i op tk+1 naar fase S, φk+1 = 2 , dan en slechts dan i als

Rik+1 > RC . ˆ Rik ≈ RC en φki = 2. De grootte van kans Pigr dat cel i over gaat naar fase G2 wordt nu bepaald door stapgrootte ∆t en de de druk die de omliggende cellen op hem uitoefenen Λki : Pigr = Pigr (∆, Λi (t)). We nemen nu een trekking uit een uniforme verdeling op [0, 1]. Voor deze stochastische variabele vi ∼ U (0, 1) geldt nu dat cel i op tk+1 naar fase S overgaat, φk+1 = 3, dan en slechts dan als: i

vi < Pig (∆t, Λi (t)). ˆ RC < Rik < Rmax en φki = 3, dan gaat deze over naar mitotische fase M, φk+1 = 4 dan i en slechts dan als:

Rik+1 > RM . ˆ Rik > Rmax en φki = 4. Dan is cel i mitotisch en geldt dus γi = 0 en ψik = 1. Deze cel gaat over tot deling als de tijd dat deze cel in de mitotische fase heeft gezeten gelijk is aan TM . Voor dochtercellen ia en ib geldt dan automatisch: k+1 φk+1 i,a = 1 en φi,b = 1.

4.2.3

Sterfte en deling

Zoals bekend uit Paragraaf 3.4.3, wordt de kans dat cel i sterft binnen tijdsinterval ∆t bepaald door de grootte van ∆t en de druk die de omliggende cellen op hem uitoefenen: Piapt (∆t, Λt ). Dit wordt op een identieke manier ge¨ımplemeerd als kans Pigr hierboven: op iedere tijdstap tk wordt een trekking van N random variabelen uit een uniforme verdeling op [0, 1] getrokken. Voor ieder van de cellen i ∈ {1, ..., N } geldt nu dat deze sterft als voor trekking ui ∼ U (0, 1) geldt: ui < Piapt (∆t, Λi (t)).

29

Numerieke implementatie Stel nu dat op tijdstip tk cel p ∈ {1, ..., N } sterft en cel q ∈ {1, ..., N }\p deelt tot cellen qa en qb . Dan geldt voor de locatiematrix Xk+1 op tijdstip tk+1 :  k+1  r1  .   ..     rk+1   p−1   k+1   rp+1   .  k+1  X =  ..  .  k+1   rq,a   .   .   .   k+1   rN  rk+1 q,b

4.3

Parameterkeuze

De parameters die voor de uiteindelijke simulatie zijn gebruikt zijn gegeven in tabel 4.1 . Deze zijn onder te verdelen in parameters die betrekking hebben op de morfologie en migratie van de cel en parameters die te maken hebben met de levensloop van de cel en de dynamiek binnen een kolonie met gezonde cellen en tumorcellen.

4.3.1

Morfologie en fysica van de cel

De morfologie-parameters zijn gebaseerd op verschillende studies. De orde van de grootte van cellen uit het long-epitheel komt uit onderzoek van Fuchs et al. [27]. De andere parameters zijn afkomstig uit [1]. Parameter R0 RC RM Es Ec cf Fb βi

Betekenis Straal van dochtercel gelijk na deling Straal cel tijdens S-fase Straal van cel vlak voor deling Stijfheid van het substraat [9] Stijfheid van de cel [9] Frictieco¨effici¨ent [23] Opwaartse trekkracht van cel [25] Mobiliteitsparameter [1]

Waarde 3.0 4.0 5.0 5.0 0.5 0.2 1.0 10

Eenheid µm µm µm kPa kPa nN s−1

Tabel 4.1: Gebruikte waarden van de celbiologische parameters.

4.3.2

Levenscyclus van de cel

De parameters die betrekking hebben op de celcyclus en de dynamiek binnen een kolonie zijn op een kwalitatieve manier gekozen. De ontwikkeling van een tumor tot een grootte van 1000 cellen zou in het echt ongeveer log2 (500) · 100 ≈ 900 dagen ≈ 8 · 107 seconden duren. Om de duur van de simulatie enigszins te verkorten is gekozen om de variabelen die ´ en seconde te maken hebben met de celcyclus en de vorming van kanker terug te schalen. E´ 30

Parameterkeuze komt hierdoor overeen met ongeveer een halve dag. Hierdoor groeit de initi¨ele tumorcel in een tijdsbestek van enkele minuten exponentieel uit tot een grote kolonie. De gebruikte waarden zijn te vinden in tabel 4.2 hieronder.

Tumor

Gezond

Type cel

Fase in celcyclus G1 S G2 M Totaal G1 S G2 M Totaal

T[s] 300 400 300 1 1001 50 0 50 1 101

Tabel 4.2: Gebruikte verblijftijden van de cellen in de diverse fasen van de celcyclus. Doordat de keuze van tijdstap ∆t = 0.1 kunnen de andere parameters nu uitgerekend worden volgens de eerder beschreven formules.

Tumor

Gezond

Type cel

Parameter T2 µ σµ Pgr (0.1) Papt (0.1) T2 µ σµ Pgr (0.1) Papt (0.1)

Betekenis Verdubbelingstijd Groeisnelheid Standaarddeviatie van groei P (φi (t + 0.1) = G2|φi (t) = S ∧ Λi (t) = 0) Drukonafhankelijke 0.1s-sterftekans Verdubbelingstijd Groeisnelheid Standaarddeviatie van groei P (φi (t + 0.1) = G2|φi (t) = S ∧ Λi (t) = 0) Drukonafhankelijke 0.1s-sterftekans

Waarde 1 · 107 2.5 5 3.33 · 10−4 6.28 · 10−5 500 20 40 1 5.15 · 10−4

Eenheid s nm/s nm/s s nm/s nm/s -

Tabel 4.3: De gebruikte groei- en deelparameters van de cellen.

31

Hoofdstuk 5

Resultaten Hieronder zullen kort enkele iteraties in het modelleerproces worden beschreven en de resultaten worden weergegeven. Veel tussenresultaten zijn in de betreffende hoofdstukken al aan de orde gekomen en zullen hier dan ook niet nog eens herhaald worden. Het uiteindelijke resultaat van dit onderzoek is te vangen in ´e´en simulatie die aan het einde beschreven zal worden.

5.1

Migratie

Voor actieve en passieve migratie is eerst gekeken naar de eenvoudige situatie met twee cellen met constante radius R = 4µm. De beginafstand tussen deze twee cellen wordt gekozen als ||r1 − r2 || = 4R = 16µm. Deze cellen botsen als ||r1 − r2 || = 2R = 8µm. Uit Sectie 3.2.2 volgt dat zonder botsing geldt M (ri ) ≈ Mi0 . Voor de cellen in dit voorbeeld geldt dan voor de situatie waarin de cellen nog niet met elkaar botsen: Vi (t) = αi M (ri ) ≈

βi R 3 0 Mi = 0.123µm/s. cf Fb2

De snelheid waarmee cel 1 en 2 elkaar naderen is dan 0.246µm/s en de cellen botsen dan ook na 8.0µm = 31.6s. (5.1) 0.2468µm/s

Als de simulatie met twee cellen wordt gedraaid, zien we dat na ongeveer 31 seconden de cellen voor het eerst met elkaar botsen. Dit komt overeen met de analytische benadering uit vergelijking (5.1). Nadat de cellen voor het eerst tegen elkaar aankomen stoten ze elkaar weer af waarna de tijdstap van ∆t = 0.1 ervoor zorgt dat de cellen tegen elkaar aan blijven botsen.

Figuur 5.1: Celmigratie.

Starttoestand eindsimulatie

5.2

Starttoestand eindsimulatie

De kolonie cellen die in deze simulatie bestudeerd wordt bestaat initieel uit 349 cellen: 347 gezonde cellen en twee centraal gelegen tumorcellen. De cellen liggen in een rooster waarbij de centra van twee naast elkaar gelegen cellen 7, 6µm van elkaar af liggen. Het cirkelvormige weefsel heeft een straal van 84µm. De straal van de cellen is uniform random verdeeld over interval [R0 , Rmax ]. De beginsituatie is weergegeven in figuur 5.2.

Figuur 5.2: Beginsituatie van simulatie; de rode punten geven de tumorcellen weer. Bij bovenstaand figuur moet worden opgemerkt dat in de grafische weergave van de simulatie de stralen van de cellen onderling identiek lijken. Dat is echter alleen een resultaat van de manier waarop de figuren zijn gegenereerd. In de simulatie zijn ze wel degelijk verschillend van grootte. Dit verschijnsel komt terug in de rest van de simulatie.

33

Resultaten

5.3

Simulatie

Er wordt nu de simulatie besproken met tijdstap ∆t = 0.1s met de parameters zoals gespecificeerd in Sectie 4.3. Vervolgens wordt een simulatie met 10.000 iteraties gedaan. Dit komt overeen met 10000 · 0.1/60 = 16, 67 minuten, oftewel 16 minuten en 40 seconden. De volgende plots zijn de resultaten na iedere minuut, dus per 600 iteraties:

t = 1 min

t = 2 min

t = 3 min

t = 4 min

t = 5 min

t = 6 min

t = 7 min

t = 8 min

t = 9 min

t = 10 min

t = 11 min

t = 12 min

t = 13 min

t = 14 min

t = 15 min

t = 16 min

Figuur 5.3: Plots van simulatie met agressieve parameterkeuze In figuur 5.3 is duidelijk te zien hoe de kolonie tumorcellen zich uitbreidt en na verloop van tijd de gezonde cellen - die aanvankelijk in aantal gelijk blijven - begint te verdringen. De parameterkeuze is dusdanig agressief dat in 16 minuten bijna het hele weefsel is overgenomen door kankercellen. Dit is geen realistisch resultaat in kwantitatieve zin. Echter zijn de parameters dusdanig gekozen dat de tijds-as wat betreft de gehele celdeling met een factor 1 dag 2 opgeschaald kan worden. Men zou dan deze simulaties kunnen interpreteren met een tijdseenheid van dagen in plaats van seconden. Bovenstaande simulatie zou dan 500 dagen vooruit gaan en in die tijd groeit tot een kolonie die ongeveer 350 maal zo groot is als oorspronkelijk. In kwalitatieve zin kan men hier dan ook goed zien hoe deze kolonie in dit tijdsframe woekert tot een gezwel. 1 seconde ∼

34

Simulatie Zoals in figuur 5.4 kan worden gezien, is hier sprake van logistische groei: in het begin neemt de populatie tumorcellen exponentieel toe, maar de groeisnelheid neemt af naarmate de kolonie dusdanig groot is geworden dat de druk van de omliggende cellen en de weefselwand ervoor zorgt dat de tumor cellen ook elkaar onderling dood gaan drukken. Uiteindelijk neemt de groei af tot een evenwicht is bereikt. Er is dus een maximum aantal cellen dat zich in het weefsel kan bevinden en dat ligt hier zo rond de 700 cellen.

Figuur 5.4: Logistisch verloop van celpopulatie in de kolonie. Verder worden de sterftekans en delingskans van de weefselcellen respectievelijk groter en kleiner door de aanwezigheid van de vele tumorcellen. Hierdoor verschuift het evenwicht van het aantal weefselcellen. Dit geeft een afname van dit aantal weefselcellen en zo is zichtbaar hoe de tumorcellen met hun “asociale” gedrag de hele kolonie overnemen. Door de toenemende druk en de korte verdubbelingstijd van tumorcellen is uiteindelijk een licht oscillerende beweging waarneembaar in het aantal kankercellen.

35

Hoofdstuk 6

Conclusie De aan dit onderzoek ten grondslag liggende hoofdvraag luidde: “Is de ontwikkeling van een tumorkolonie realistisch te modelleren op celniveau?” Om tot zo’n model te komen moest allereerst rekening worden gehouden met de fysische eigenschappen van gezonde cellen en kankercellen, de manier waarop cellen met elkaar communiceren en de wijze waarop passieve en actieve migratie plaatsvinden. Tevens moest de gehele celcyclus inclusief groei, deling en sterfte meegenomen worden, en moest worden gekeken naar de invloed hierop door de druk die omliggende cellen uitoefenen. Ook hierin zijn de verschillen tussen gezonde cellen en kankercellen meegenomen. Tot slot zijn de parameters voor deling en sterfte statistisch en biologisch zo realistisch mogelijk beargumenteerd. Om de gehele ontwikkeling van cel tot gezwel binnen een acceptabel tijdsbestek te simuleren, was voor de delingsparameters nog wel een opschaling in de tijd nodig. Bovenstaande heeft uiteindelijk geleid tot een model waarin ´e´en gemuteerde cel uitgroeit tot een longtumor die gezonde cellen van het weefsel dooddrukt. Het aantal tumorcellen dat in de simulatie is ontstaan is niet groot genoeg om opgemerkt te worden door het organisme waarin deze tumor is ontstaan. De tumor van uiteindelijk 700 cellen- samen nog geen 0.1 millimeter in doorsnede- staat niet in verhouding tot de 50.000 cellen die nodig zijn voor een tumor die tot uitzaaiingen kan leiden. De uitkomsten van deze simulaties zijn dus niet van dien aard dat op deze schaal moet worden gekeken, maar eerder naar de beginstadia van de ontwikkeling van een tumor. Deze lijken op een realistische manier gesimuleerd te zijn. Tijdens de implementatie is de nadruk gelegd op tijdseffici¨entie, waardoor uitbreiding naar grotere celkoloni¨en en complexere modellen eenvoudig te maken zal zijn. Door de opschaling van de delingsparameters zijn de uitkomsten vooral te interpreteren in kwalitatieve zin. De resultaten zijn veelbelovend en lenen zich tot uitbreiding en verder onderzoek. Hieruit volgt een positief antwoord op de hoofdvraag: de ontwikkeling van een tumorkolonie is realistisch te modelleren op celniveau.

Hoofdstuk 7

Discussie In het model zijn een aantal aannames gedaan die wellicht bij uitbreiding van dit model nog eens goed heroverwogen moeten worden. ˆ Tijdschaal. Ten eerste is de tijdsschaal van migratie realistisch genomen maar die van de celcyclus niet. De duur van de celcyclus van een gezonde cel ligt rond de 24 uur, in plaats van de 20 minuten die in de uiteindelijk simulatie zijn genomen. De cyclus van een kankercel is weliswaar een stuk korter van duur dan die van een gezonde cel, maar zeker niet de 1,5 minuut waarvoor uiteindelijk gekozen is. Snellere implementatie (wellicht in een andere programmeertaal), langere simulatie en meer geduld zouden hierin realistischere uitkomsten kunnen geven, al moet de kwalitatieve waarde van de hier gedane simulatie zeker niet onderschat worden. ˆ Indeling celcyclus. Uit onderzoek [28] blijkt dat de verhouding van de tijdsduur van verschillende fases uit de celcyclus van longepitheel anders is dan die gebruikt is in de simulatie. De tijdsduur van G2 zou een factor acht lager horen te liggen dan die van G1. Een aanpassing hiervoor zou gemakkelijk gemaakt moeten worden. ˆ Morfologie en fysiologie tumorcellen. In dit model zijn een aantal parameters gelijk verondersteld voor tumorcellen en gezonde cellen. Dit is niet overal terecht. Zo groeien tumorcellen niet alleen harder en ongecontroleerder, maar zijn ze ook typisch groter dan gezonde cellen. Ook verschillen ze in elasticiteit en dus in deformatie en bewegingssnelheid. Onderzoek naar betere kwantificering van deze parameters zou gewenst zijn. ˆ Drukafhankelijkheid. In dit model wordt de kans op groei en de kans op deling van een cel alleen be¨ınvloed door de direct aangelegen cellen, terwijl het wellicht realistischer is als alle nabijgelegen cellen meegenomen worden. Immers, een cel ontvangt ook hiervan de signalen en reageert hierop. ˆ Stabiliteit. In dit onderzoek is geen theoretische onderbouwing van de stabiliteit gedaan. Trial-and-error gaf uiteindelijk bevredigende resultaten, maar vanuit wiskundig oogpunt zou een goede onderbouwing hiervan zeker niet misstaan. ˆ Nauwkeurigheid. De nauwkeurigheid van de Euler Voorwaartse methode is slechts
O ∆t . Een hogere orde methode zou de tijdsintegratie van de celgroei en migratie nauwkeuriger maken.

Discussie

7.1

Aanbevelingen

Uit de conclusie komt naar voren dat het model zoals dat hier is gepresenteerd en ge¨ımplementeerd goed bruikbaar is voor uitbreiding. Vanuit de discussie zijn een aantal punten naar voren gekomen die wellicht aangepast zouden moeten worden om een nog realistischer model te hebben. Hieronder worden een aantal aanbevelingen gedaan die wellicht kunnen aanzetten tot en bijdragen aan een nuttig vervolg van dit onderzoek. Hierbij is het zinvol samenwerking aan te gaan met een faculteit of onderzoeksinstituut in de celbiologie of oncologie. ˆ Analyse stabiliteit en nauwkeurigheid. Een grondige analyse van de stabiliteit meer vertrouwen geven in uitkomsten van gedane simulaties en wellicht kunnen leiden tot versnelde simulatie door een grotere stapgrootte ∆t. ˆ Aanpassing parameters. Een uitgebreidere biologische onderbouwing van de parameters, eventueel in samenwerking met cel-biologen en kankeronderzoek ter verkrijging van een realistischer model. ˆ Koloniegrootte. Bij een longtumor vinden uitzaaiingen plaats bij tumorgrootten in de orde van 50.000 cellen. Een kolonie van enkele honderden cellen is dan ook een erg kleine schaal. Wellicht zou met meer rekenkracht en -tijd, en wat aanpassingen in de implementatie, dit model met een factor 100 opgeschaald kunnen worden, waardoor uitzaaiingen gemodelleerd kunnen worden. Voor nog grotere tumoren is het aan te raden over te stappen naar continue modellen. ˆ Zuurstofconcentratie. De groei van cellen wordt niet alleen be¨ınvloed door omliggende cellen maar ook door fysische factoren zoals het zuurstofgehalte. Hier zijn al verscheidene onderzoeken naar gedaan, mede door medestudente Liselot Arkesteijn. Wellicht is het mogelijk deze onderzoeken te combineren tot een model voor de ontwikkeling van een tumor waarbij ook rekening wordt gehouden met de sterfte van kankercellen door zuurstofgebrek. ˆ Angiogenese. Binnen een tumor vindt angiogenese plaats, een proces waarbij bloedvaten worden gevormd om meer zuurstof naar de tumorcellen toe te voeren. Dit zou in samenwerking met bovenstaand onderzoek naar de invloed van de zuurstofconcentratie kunnen worden meegenomen in het model. ˆ Uitzaaiing. Een uitzaaiing vindt plaats als een tumor dusdanig gegroeid is, dat de druk van omliggende cellen de tumorcellen door de wand van het weefsel heen drukken en zo op een andere plaats in het lichaam terechtkomen. Als er naast het weefsel een bloedvat is, kan de tumorcel ergens anders in het lichaam via transport door de bloedbaan een uitzaaiing veroorzaken. Enkele eenvoudige aanpassingen zouden voldoende moeten zijn om een bloedvat te construeren en een tumorcel hier naar binnen te laten drukken bij genoeg druk Λi van omliggende cellen. ˆ Samendrukken nabijgelegen orgaan. Een gevaar bij tumoren is dat ze dusdanig in grootte toe kunnen nemen dat nabijgelegen organen worden verdrukt en in hun functie worden belemmerd. Als een orgaan naast het weefsel in de simulatie zou worden geconstrueerd zou dit op celniveau een model geven van dit proces. ˆ Weefselkeuze. Door aanpassing van parameters zou eenvoudig de keuze van het weefsel kunnen veranderen van longepitheel naar een vergelijkbaar weefsel. ˆ 3D-model. Een interessante uitbreiding zou het omschrijven naar een 3D-model zijn. Dit is grafisch wel een stuk uitdagender, maar wat implementatie betreft is dit relatief eenvoudig.

38

Literatuurlijst [1] F.J Vermolen, A.Gefen: A semi-stochastic cell-based formalism to model the dynamics of migration of cells in colonies,Biomechanics Modeling in Mechanobiology, 11(102): 183-195 (2012). [2] Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, Peter Walter: Molecular Biology of the Cell, Fourth Edition, New York: Garland Science (2002). [3] Vickaryous MK, Hall BK: Human cell type diversity, evolution, development, and classification with special reference to cells derived from the neural crest, Biological reviews, Volume: 81, Pages: 425-455 (2006). [4] http://cscbiologie.jouwweb.nl, bezocht op 7 januari 2013. [5] Alberts et al.: Essential Cell Biology, Second Edition, 2003. [6] http://www.biology-online.org/, bezocht op 7 januari 2013. [7] Ananthakrishnan R, Ehrlicher A: The Forces Behind Cell Movement, International Journal Biological Science; 3(5):303-317 (2007). [8] Abercrombie M. The Croonian Lecture, The Crawling Movement of Metazoan Cells. Proc Roy Soc Lond B. 1980;207:129147 (1978). [9] Reinhart-King, C.A., Dembo, M., Hammer D.A.: Cell-cell mechanical communication through compliant substrates, Biophys. J., 95, 60446051, (2008). [10] Aleksander S. Popel: Biomechanical measurement in cells, Department of Biomedical Engineering, Whitaker Biomedical Engineering Institute, School of Medicine, Johns Hopkins University, Baltimore, Maryland [11] Anand P., Kunnumakkara A.B., Kunnumakara A.B., Sundaram C., Harikumar K.B., Tharakan S.T., Lai O.S., Sung B., Aggarwal B.B. Cancer is a Preventable Disease that Requires Major Lifestyle Changes. Pharm Res. 2008 September; 25(9): 20972116. [12] Alberg AJ, Samet JM (2010). ”46”. Murray & Nadel’s Textbook of Respiratory Medicine (5th ed.). Saunders Elsevier. ISBN 978-1-4160-4710-0. [13] Ferlay J, Shin HR, Bray F, Forman D, Mathers C and Parkin DM. GLOBOCAN 2008 v1.2, Cancer Incidence and Mortality Worldwide: IARC CancerBase No. 10 [Internet]. Lyon, France: International Agency for Research on Cancer, (2010). [14] Jemal A, Bray, F, Center, MM, Ferlay, J, Ward, E, Forman, D (February 2011): ”Global cancer statistics”. CA: a cancer journal for clinicians 61 (2): 6990. [15] Sasco AJ, Secretan MB, Straif K. :Tobacco smoking and cancer: a brief review of recent epidemiological evidence. Lung Cancer 45 (Suppl 2): S3-9. PMID 15552776, (2004).

Literatuurlijst [16] Doll R, Hill AB. A study of the aetiology of carcinoma of the lung. BMJ 1952;ii:127186. [17] Charles J. Sherr, Cancer Cell Cycles, SCIENCE, VOL. 274, 6 december (1996) [18] Subra Suresh, Biomechanics and biophysics of cancer cells. 02-04-2007 [19] Igawa S et al.: Nitric oxide generated by iNOS reduces deformability of Lewis lung carcinoma cells. Cancer Sci;95:3427, (2004). [20] Werner R. Loewenstein and Yoshinobu Kanno: Intercellulair communication and tissue growth, J Cell Biol 33:225-234. Published May 1, 1967, doi:10.1083/jcb.33.2.225, (1967). [21] Gefen, A, Department of Biomedicine of Tel Aviv University. [22] http://mens-en-gezondheid.infonu.nl/ziekten/24856-hoe-snel-groeien-kankercellen.html, bezocht op 29 januari 2013. [23] Gefen, A.: Effects of virus size and cell stifness on forces, work and pressures driving membrane invagination in a receptor-mediated endocytosis, Journal of Biomechanical Engineering (ASME), 132, 084501-1-084501-5-2010. [24] J. Michael Steele. Stochastic Calculus and Financial Applications, Springer-Verlag New York, (2001). [25] Lemmon, C.A., Chen, C.S., Romer, L.H.: Cell traction forces direct fbronectin matrix assembly, Biophys. J., 96, 729-738, (2009). [26] Wang, JHC., Lin, J.-s.: Cell traction force and measurement methods, Biomech. Model Mechanobiol., 6, 361-371, (2007). [27] Fuchs et al.: Differentiation of human alveolar epithelial cells in primary culture: morphological characterization and synthesis of caveolin-1 and surfactant protein-C., Cell and Tissue Research, 2003 Jan; 311(1):31-45. [28] Bruce et al.: Cell size, cell cycle, and α-smooth muscle actin expression by primary human lung fibroblasts, Lung Physiol November 1, 1998 vol. 275 no. 5 L998-L1005, (1998).

II

Appendix A

Appendix A.1 A.1.1

%% % % % % % % % % % % % % % % % % % % % %

Implementatie in matlab Simulatie

Implementation of 'A semi−stochastic cell−based formalism to model tumorgrowth'. %% Date : 29−12−2012 Institution: TU Delft, Faculty of Electrical Engineering, Mathematics & Comuter science Part of : BSc Thesis Delft institute of Applied Mathematics By : Rune van der Meijden Supervisor : Dr. Ir. F.J.Vermolen Sketch of the situation: Model of epithelial tissue from the lungs. The tissue is circular with radius 'RadColony' and has a boundary that has the same elastic properties as the epithelium itself. In the middle of the tissue, a healthy cell turns into a benign tumorcell and a lung carcinoma starts growing. Cell morphology: Cells are assumed to be hemi−spherical with radius R0 <= R(i) <= Rmax. Location of each cell is given by it's middlepoint r(i). Startcondition: The colony exists of (n1=)347 healthy cells in a circular raster and (n2=)2 benign tumor cell in the middle.

clear all; close all;

%% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % −−>>> Parameters <<<−− % −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Integration dt=0.25; % Time [s] % \\

Cellmorphology // R0=3e−6; Rc=4e−6; Es=5e3; Ec=0.5e3; Fh=1e−9; F=Fh; Cf=0.2; beta=10;

Rm=5e−6;

% Radii [m] % Elasticity [kPa] % Traction force[N] % Mobility Parameters

% \\ Cell cycle phases health( H) and tumor ( C)−cells // TG1 H = 300; TG1 C= 50; % Duration of Gap 1 (=Gap 2) TS H = 250; TS C = 0; % Duration of Synthesis

Appendix TMit = 1; TCell H = (2* TG1 TCell C = (2* TG1 T2 H = TCell H * T2 C = TCell C * % \\

% \\

% \\

H + TS H + TMit); C + TS C + TMit); 400; 1.5;

Growth // GrowSpeed H = (Rm−R0)/(2* TG1 H); GrowSpeed C = (Rm−R0)/(2* TG1 C); sigma H = GrowSpeed H*2; sigma C = GrowSpeed C*2; Pgrow H = dt/TS H; Pgrow C = dt/TS C; Death // Papt H = 1−2ˆ((1/T2 H − 1/TCell H)*dt); Papt C = (1−2ˆ((1/T2 C − 1/TCell C)*dt)); Colony // ColonyRadius = 21*Rc;

% Duration of Mitosis % Duration of total cell−cycle % Doubling time of colony

% Average speed of a growing cell % Standard deviation of growth % Prob. S−>G2 in dt, no pressure

% Prob. of apoptose, no pressure

% Radius of colony

% \\ Pressure function // fPress = @(var) 10./(1 + 9*exp(var)) ;

%% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % −−>>> Start Conditions <<<−− % −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % \\

Location // x h = importdata('xyVec.mat'); x c = [0,0;x h(174,:)]; x h(174,:) = []; x=[x h;x c]; [n1,m]=size(x h); [n2,m]=size(x c); N=n1+n2;

% \\ Size // %load StartFase.mat R=[R0 + (Rm−R0)*rand(n1,1);Rc;Rm]; TumorIndicator=[zeros(n1,1);ones(n2,1)]; % \\ Movability // M0=Fhˆ2/(2*piˆ2*Es*R0ˆ4); eps=M0*exp(−Es/Ec*30e−6/R0); a=beta*R0ˆ3*F/(Cf*Fhˆ2); lbd=Es/Ec;

% Location of healthy cells % Location of tumor cells % Location of all cells % Number of healthy(n1) and tumor(n2) cells % Total number of cells

% Phases of cells % Random size of cells % Tumor Indicator: 0 = healthy, Tumor = 1;

% % % %

M jˆ0 assumed constant Detection Treshold alpha lambda

% \\ Phase of Cell Cycle : Gap 1 = 0, Synthesis= 1, Gap 2=2, Mitosis=3 // CellCyclephase = 2*(R > Rc & R < Rm); Mitosis = zeros(N,1); % \\ Plot & Movie stuff // CellSizePlot = 10; as=32*[−R0 R0 −R0 R0]; MaxSimulationLength = 10000;

% Size of cells in plot (constant) % Plot size % [Frames];

% \\ Cell administration // HCells=nan(MaxSimulationLength,1); CCells=HCells; TCells=HCells+CCells; MaxLambda = HCells; MeanLambda = HCells; NumLambda = HCells;

%% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

IV

Implementatie in matlab % −−>>> Simulation <<<−− % −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− Index=1; h = waitbar(0,'Duurt laaang...')

while (n2 <= N) & (Index < MaxSimulationLength) %i = 1 : length(t) waitbar(Index/MaxSimulationLength,h) % \\ Plot // figure(1); plot(x((TumorIndicator==0),1),x((TumorIndicator==0),2),'o','MarkerSize',... CellSizePlot); hold on;plot(x((TumorIndicator==1),1),x((TumorIndicator==1),2),'o',... 'MarkerFaceColor','r','MarkerSize',CellSizePlot); rectangle('position',[−ColonyRadius −ColonyRadius 2*ColonyRadius 2*ColonyRadius],'curvature', hold off; axis(as);

% \\ Record movie frame every minute // frames(Index) = getframe; if mod(Index,60/dt) == 1 %Figurename = ['C:\Users\Jorrit\Desktop\Rune Bep\Plaatjes\', 'SituatiePlot na ', num2str( Figurename = ['C:\Users\Rune\Dropbox\BachelorProject\Plaatjes\', 'SituatiePlot na ', ... num2str(dt*(Index−1)/60), ' minuten.jpg']; saveas(figure(1), Figurename) end % \\ Cell pressure // Dist = DistMat(x); % Dist(i,j)= | | rj − ri | | RadSum = CumulativeRadius(R); % RadSum(i,j) = Ri + Rj hMat = max(RadSum−Dist,0)/2; % hMat(i,j) is impinging.. hMat(find(eye(N)))=zeros(N,1); % ..of cell i because of cell j hWall = max(0,NormOfVec(x)+ R − ColonyRadius); % Impinging from tissuewall MPress = 4*Ec/(15*sqrt(2)*pi)*(2*hMat./RadSum).ˆ(5/2); % Pressure from cells MPressWall = 4*Ec/(15*sqrt(2)*pi)*(hWall./R).ˆ(5/2)'; % Press from tissuewall PressureSum = sum(MPress)+ MPressWall;

% \\ Growing // fLambda = 1.2*fPress(1*PressureSum)'; % f(Lambda i), see thesis Pgrow = Pgrow H*(1−TumorIndicator).*fLambda + Pgrow C*TumorIndicator; % Grow prob.(press) Growing =(CellCyclephase==0)+(CellCyclephase==2); % Only growth in G1 & G2 randVec = randn(N,1); R = R + dt*((1−TumorIndicator).*(Growing.*GrowSpeed H+randVec* sigma H/sqrt(dt))... + TumorIndicator.*(Growing.*GrowSpeed C+ randVec* sigma C/sqrt(dt))); % \\ Cell Cycle phase changes // CellCyclephase = CellCyclephase+(min((CellCyclephase==0),(R > Rc))); %G1 −> Synth CellCyclephase = CellCyclephase+(min((CellCyclephase==1),(rand(N,1)<=Pgrow))); %Synth−> G2 CellCyclephase = CellCyclephase+(min((CellCyclephase==2),(R>Rm))); %G2 −> Mit % \\ Apoptosis // Psurvive = (1−Papt H).ˆ(1/1.2*(fPress(−1*PressureSum))).*(1−TumorIndicator')... +(1−Papt C).ˆ(1/1.2*(fPress(−.05*PressureSum))).*TumorIndicator'; % Apoptosis prob with press Papt = (1 − Psurvive'); if (n2==1); Papt(logical(TumorIndicator))=0; end %(Optional) Avoiding n2 = 0 DeathIndex = find(rand(N,1) <= Papt); % Index for dying cells x(DeathIndex,:)=[]; R(DeathIndex)=[]; TumorIndicator(DeathIndex)=[]; CellCyclephase(DeathIndex)=[]; Mitosis(DeathIndex)=[]; n1=sum(1−TumorIndicator); n2=sum(TumorIndicator); N=n1+n2; % \\ Mitosis & Migrating? //

V

Appendix Mitosis = Mitosis+(CellCyclephase==3)*dt; MoveIndex = find(Mitosis < TMit); DivIndex = find(Mitosis >= TMit); % Dividing cells MigrationIndex = min((TumorIndicator==0),(Mitosis==0));

% Time cell is in Mitosis phase % Passive movement % Active movement

% \\ Strain energy density // Dist = DistMat(x); RadSum = CumulativeRadius(R); MjRi = exp(−lbd/Rc*Dist)*diag(M0); % Exponential decay of signal MjRi = MjRi.*(MjRi>=eps); % Epsilon := detection treshold Mstrain= MjRi*diag(MigrationIndex); % Only non−mitosis non−tumor cells 'feel' strain hMat=max(RadSum−Dist,0)/2; hWall= max(0,NormOfVec(x)+R−ColonyRadius); MPress = 4*Ec/(15*sqrt(2)*pi)*(2*hMat./RadSum).ˆ(5/2);% Pressure from cells MPressWall = 4*Ec/(15*sqrt(2)*pi)*(hWall./R).ˆ(5/2); % Pressure from tissuewall M = Mstrain−MPress; M(find(eye(N))) = zeros(N,1); Msum = M0 + sum(abs(M))' + MPressWall; % \\ Direction of movement // [Vx Vy] = DirectionVec(x); % Vector from cell i to cell j zx=diag(M*Vx'); zy=diag(M*Vy'); % Decompose in x and y direction zWall = −x.*repmat(MPressWall./NormOfVec(x),1,2); % if (sum(any(x,2))˜=N); zWall(˜any(x,2),:) = 0; end z=[zx, zy] + zWall; % Direction of movement zh= z ./ [NormOfVec(z) NormOfVec(z)]; % DirectionVec(z); if (sum(any(z,2))˜=N); zWall(˜any(z,2),:) = 0; end % \\ Movement // if length(MoveIndex) > 0 x(MoveIndex,:)=x(MoveIndex,:)+dt*a*[Msum(MoveIndex),Msum(MoveIndex)].*zh(MoveIndex,:); end

% \\ Mitosis // vDiv = RandUnitVec(length(DivIndex)); xNew = x(DivIndex,:)+ dt*a*[Msum(DivIndex),Msum(DivIndex)].*(zh(DivIndex,:))−.8*R0*vDiv; x(DivIndex,:)=x(DivIndex,:)+ dt*a*[Msum(DivIndex),Msum(DivIndex)].*(zh(DivIndex,:))+0.8*R0*vD x = [x ; xNew]; R(DivIndex) = R0; R = [R ; R(DivIndex)]; CellCyclephase(DivIndex) = 0; CellCyclephase = [CellCyclephase ; CellCyclephase(DivIndex)]; Mitosis(DivIndex) = 0; Mitosis = [Mitosis; Mitosis(DivIndex)]; TumorIndicator = [TumorIndicator ; TumorIndicator(DivIndex)]; n1 = sum(1−TumorIndicator); n2 = sum(TumorIndicator); N = n1 + n2; % \\ Cell administration // HCells(Index)=n1; CCells(Index)=n2; TCells(Index)=N; TumorIndex = (TumorIndicator == 1); Index = Index + 1; display(N); end

%% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % −−>>> Plotting <<<−− % −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % \\ Plot cell population // figure(2); plot((1:Index)*dt,HCells(1:Index));hold on; plot((1:Index)*dt,CCells(1:Index),'r'); plot((1:Index)*dt,TCells(1:Index),'k'); xlabel('Tijd [s]'); ylabel('Cel populatie in kolonie'); legend('Gezonde cellen', 'Tumor cellen', 'Totaal aantal cellen') hold off;

VI

Implementatie in matlab

%% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % −−>>> Make Movie <<<−− % −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− Datenow = date; Timenow = clock; Timenow = Timenow([4 5]); PathName = 'C:\Users\Rune\Desktop\BEPfilmpjes\Simulation of tumor growth '; TimeStamp = [Datenow , ' ' num2str(Timenow(1)) ,'', num2str(Timenow(2))]; VideoNameString1 = [PathName, TimeStamp,'.avi']; writerObj1 = VideoWriter(VideoNameString1); writerObj1.Quality = 55; writerObj1.FrameRate = 60; whitebg('w'); open(writerObj1); open(writerObj1); writeVideo(writerObj1,frames); close(writerObj1); VideoNameString2 = [PathName, 'resolution 1to10) ', TimeStamp,'.avi']; writerObj2 = VideoWriter(VideoNameString2); writerObj2.Quality = 55; writerObj2.FrameRate = 20; whitebg('w'); open(writerObj2); open(writerObj2); writeVideo(writerObj2,frames(1:20:end)); close(writerObj2);

VII

Appendix

A.1.2

Subfuncties

function[ Y ] = DistMat( X ) %Returns matrix with distance between points % Detailed explanation goes here R=X; %Matrix X has [x1 y1; x2 y2; ..] with points in R2: r1=(x1,y1) [m,n]=size(R); Ax=ones(m,1)*R(:,1)'; Bx=(Ax−Ax'); Cx=Bx.*Bx; Ay=ones(m,1)*R(:,2)'; By=(Ay−Ay'); Cy=By.*By; D=sqrt(Cx+Cy); Y=D;

end

function [CumRad] = CumulativeRadius(Rvec) % Input: vector with Radii of m cells % Output: Matrix CumRad with CumRad(i,j)=Rvec(i)+Rvec(j) R=Rvec; [m,n]=size(R); RR=repmat(R,n,m); CumRad=RR+RR'; end

function [ NormVec ] =NormOfVec( x ) %NormOfVec returns a vector of 2−norms of the rows of x % x = [x1 y1; .. .. xn yn] −> NormOfVec(x) = [ | | ( x1,y1 ) | | ; ...; | | ( xn,yn ) | | NormVec = sqrt(sum(x.ˆ2,2)); end

VIII

Implementatie in matlab function [Vx Vy] = DirectionVec(r ) %DirectionVec returns normalvector V {ij} from eq (11) %% % The input is an array r=[r1;r2;..;rm]=[rx,ry] with ri=[xi,yi], a point in % R2−space. This programm returns 2 matrices Vx and Vy with % Vx(i,j) is the x−coordinate of normalvector ri −> rj % Vy(i,j) is the y−coordinate of normalvector ri −> rj R=r; [m,n]=size(R); %Compose matrix Dx= Dx(i,j)=xj−xi X=ones(m,1)*R(:,1)'; %[rx ; rx;...;rx] Dx=(X−X'); %Dx(i,j)=xj−xi Y=ones(m,1)*R(:,2)'; %Idem for y.. Dy=(Y−Y'); %Use Distmat.m to compes D(i,j ) = | | rj−ri | | DM=DistMat(R); %Compose V Vx = Dx./DM; Vy = Dy./DM; Vx(isnan(Vx)) = 0; Vy(isnan(Vy)) = 0;

%Deviding by 0 means NO normalvector

end

function [ v ] = RandUnitVec( n) %Returns unitvector in uniform random direction theta=2*pi*rand(n,1); v=[cos(theta), sin(theta)]; end

IX