Table Of Content

Statistical mechanics model of angiogenic tumor growth Antońio Luis Ferreira,1 Dorota Lipowska,2 and Adam Lipowski3 1Departamento de Fisica and I3N, Universidade de Aveiro, 3810-193 Aveiro, Portugal 2Faculty of Modern Languages and Literature, Adam Mickiewicz University, Poznan´, Poland 3Faculty of Physics, Adam Mickiewicz University, Poznan´, Poland We examine a lattice model of tumor growth where survival of tumor cells depends on the supplied nutrients. When such a supply is random, the extinction of tumors belongs to the directed 2 percolation universality class. However, when the supply is correlated with distribution of tumor 1 cells, which as we suggest might mimick theangiogenic growth, the extinction shows different, and 0 most likely novel critical behaviour. Such a correlation affects also the morphology of the growing 2 tumorsand drastically raise tumorsurvival probability. n a PACSnumbers: 87.18.Hf J Keywords: tumorgrowth,angiogenesis, absorbingstates 0 1 Duetoseriesofmutationsofgenesresponsibleforpro- hopethatfurtherdevelopmentwillresultinlatticemod- ] liferation,somecellsmightinitiatetheabnormalandun- els taking into account some other important, but so far h controlled growth process, commonly named as tumor neglected, factors like heterogeneity,immune cells or the c e [1]. In the first stage of the process, called avascular role of chemoattractants. m growth,thereisnobloodsupplyandthegrowthislimited Discrete lattice modeling allows us to describe tumor - by the amount of oxygen and nutrients that the tumor growth problem using tools developed in statistical me- t a receivesthroughitssurface. Thegrowingtumorinitiates, chanics for studying complex systems. A distinctive fea- st however,anumberofaccompanyingprocessesinthehost ture of these systems is spontaneous emergence of cer- t. environment. In particular, large demand for nutrients tain properties, which cannot be traced to the character a stimulatestumorcellstoproduceangiogenicfactorsthat of individual parts. It is believed that life, conscious- m regulate the formation and growth of new blood vessels ness, or functioning of ant colonies are examples of such - in the region. This stage is called angiogenesisand it in- emergent phenomena [13]. There is an increasing evi- d n termediatesbetweenavascularandvasculargrowth. The dence that cancer can also be considered as an emergent o vascular growth is the third stage, which begins when property, and thus, developing statistical-mechanics ap- c the blood vessels have reached the tumor. In this stage proaches seems to be very promising [14–16]. In partic- [ the tumor receives a vast amount of nutrients and can ular, one can hope that such bottom-up modeling will 2 grow much larger than it was possible during the avas- helpus toexplaintumorgrowthintermsofcellparame- v cular growth. Moreover, the vasculature might be used ters, which might contribute to its better prediction and 0 to spread tumor cells throughout the body of the host, control. 3 which very often leads to its death. 5 The ultimate goal, namely builiding realistic and 3 Tumor growth is a very complex process and to fully testable against real data models, will most likely re- 0. understand its nature, one has to resort to computa- quire very complex, multi-scale models, which will be 1 tionaltechniques,whichwouldsupplementbiologicaland difficult to understand without extensive computer sim- 1 medical approaches [2–6]. Various models were used to ulations. To develop some intuitive understanding of tu- 1 describe tumor growth [7]. Initially, they were contin- morgrowth,itis thus desirableto examinesome simpler : v uous models formulated in terms of partial differential models, which hopefully contain important ingredients i equations and studied mainly from mathematical point of the process. In the present Letter, we examine a sim- X of view [8]. More recently, inclusion of biomechanical ple lattice growth model where tumor cells survival and r a details and coupling of tumor growth with the vascular- breeding depend on the supplied nutrients. Mimicking ization process shifted modeling toward more physical angiogeneticprocesses,weassumethatthe supplyofnu- approach [9]. trients is positively correlated with the distribution of An alternative to differential-equation approach is tumorcells. Itturnsoutthatsuchacorrelationsubstan- based on discrete lattice models such as, for example, tially changes the statistical mechanics behaviour of the nonequilibriumQ-statePottsmodels. Withsucha mod- model. In particular, the extinction of tumor does not eling,onecanimplementseveralaspectsofcelldynamics, belong to the expected directed-percolation universality whicharedifficultto treatsimultaneouslyusingcontinu- classandthecriticalexponentβ describingtheorderpa- ous modeling, as e.g. multiplication, competition, aging, rametermostlikelygetstheclassicalvalueβ =1,evenin death, mutations and even adhesion or chemotaxis [10]. thed=1versionofthemodel. Suchacorrelationaffects Somemodelsofthiskindweresuccessfullytestedagainst also the morphology of growing tumors and drastically clinical data of certain forms of cancer [11, 12]. One can raises the tumor-survival probability. 2 Inourmodel,eachsiteofad-dimensionallatticeeither 1 is occupied by tumor, nutrient, both tumor and nutri- 0.9 ent, or is empty. At a rate p, nutrients are supplied to 0.8 a chosen site of the lattice, provided that the site is not occupiedalreadybyanutrient. Theroulette-wheelselec- 0.7 tion [17] is used to choose the site for such a supply and 0.6 tiopshfaeorianccmccourreeprtaeieessrdpedo∆bnynd>uitnturg0miewtnoaetrksig(ewhssutip=nwptol1dye+apdce∆ucnoe)dustonortonnfaoonwrtgmhi(oeawgttheio=nenric1to)he.feffTnesceihttwees x,xtn 000...345 xxxxnntt,,,, ∆∆∆∆====0303 g(x)10t--10-- ..21055 ββ==1 0.0.276 lo-2.5 blood vessels in the vicinity of tumor cells. At a rate 0.2 -3 1−p,a tumorcellonarandomlychosensite is updated. -3.5 0.1 -4 -3.5 -3 -2.5 -2 -1.5 -1 The tumor cell survives provided that there is a nutri- log (p-p) 0 10 c ent on this site, otherwise it dies. The surviving tumor 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 cell consumes the nutrient and attempts to breed pro- p vided that there is a site without a tumor cell among its FIG.1. Steady-statedensitiesofnutrients(xn)andoftumor nearest neighbours. Closely related models but without cells (xt) as a function of p calculated for d = 1 and ∆ = 0 any preference in nutrient supply (∆ = 0) were already and3. Theinsetshowsthebehaviourofxtclosetothecritical point(log-logscale). Whilefor∆=0thedirected-percolation studied [18]. scaling is seen, muchdifferent behaviour appears for ∆=3. To examine our model we used Monte Carlo simulations[19]. Forvariouspand∆,wemeasuredthesteady- statedensitiesoftumorcellsxt andofnutrientsxn. Sim- -0.4 ulations started from a random initial configuration and -0.6 the model relaxed until a steady state was reached. We examined lattices of various sizes N to ensure that the -0.8 obtained results are N-independent. We also measured -1 tthheetuimnietdoefpteinmdeenccoerroefstphoendtuinmgortoceNllduepndsiattyexatt(tt)e,mwpittsh. [x(t)]10t -1.2 DP The density xt(t) for each t is an average over indepen- og l -1.4 dent runs. At the critical point, the density xt(t) is expected to havea power-lawdecayxt(t)∼t−δ, whereδ is -1.6 a characteristic exponent [20]. -1.8 First, we describe results for the d = 1 version of the model. Simulations show that, for sufficiently large p, -2 1 2 3 4 5 6 the modelremains in anactive phase with xt >0,which log (t) 10 terminates at a critical point pc depending on ∆. For FIG. 2. Time dependence of the density of tumor cells xt(t) p < pc, the steady state of the model is an absorb- (log-log scale) for d = 1, ∆ = 0 and (from top) p = 0.37, ing state xt = 0 and xn = 1 (tumor cells die out of 0.3695, 0.3692, 0.36905, 0.3689 0.368, 0.366, and 0.36. The the lack of nutrients). At p = pc, the model undergoes dotted straight line has the slope corresponding to the di- the phase transitionfrom an active into absorbingphase rected percolation valueδ=0.159. and we expect that xt is the corresponding order parameter. As it is already known, models with a single absorbing state are expected to belong to the so-called is obtained from the behaviour of xt(t) (Fig. 3). How- directed-percolation (DP) universality class [20]. In this ever, at p = pc the decay seems to be described by the universality class and for d = 1, the critical exponent exponentδ =0.60(5),whichisagainmuchdifferentthan δDP = 0.159(1) and the decay of the order parameter theDPvalue0.159. Incalculationsfor∆=0 and 3,and upon approaching the critical point is described by the close to critical points, systems as large as N = 5x105 exponent βDP = 0.276(1). The calculated steady-state were used and simulation and relaxation times were of values of xt and xn for ∆ = 0 and 3 are presented in the order of 106. Further away from critical points, less Fig.1. For ∆ = 0, we estimate pc = 0.3691(2) and the extensive simulations were required. exponentβ isveryclosetothe DP value (insetinFig.1). For ∆=0.1 on a shorter time scale (t∼102), one can The decay of xt(t) at the critical point (Fig.2) confirms notice a slower, DP-like decay of xt(t), that, however, the DP universality class in this case. turnstomuchfasterdecayonalongertimescale(Fig.4). Much different behaviour is seen for ∆ = 3. In this Such a behaviour indicates the proximity of the DP be- case, we estimate pc = 0.1765(2) and the exponent haviour that occurs at ∆ = 0. However, critical expo- β = 1.0(1) (inset in Fig.1). The same estimate of pc nents β and δ for ∆ = 0.1 are nearly the same as for 3 -0.5 0 -1 -0.5 -1.5 DP -1 [x(t)]0t -2-.25 [x(t)]0t -1.5 0 .02.53 og1 og1 xt 0.2 l -3 l 0.15 -2 δ=0.6 0.1 -3.5 δ=0.6 0.05 -2.5 -4 0 0.355 0.36 0.365 0.37 p -4.5 -3 1 2 3 4 5 6 0 1 2 3 4 5 6 7 log (t) log (t) 10 10 FIG. 3. Time dependence of the density of tumor cells xt(t) FIG. 4. Time dependence of the density of tumor cells xt(t) (log-logscale)ford=1, ∆=3and(fromtop)p=0.19,0.18, (log-log scale) for d =1, ∆ =0.1 and (from top) p= 0.356, 0.178, 0.177, 0.1765, 0.176, 0.175, 0.17. The dotted straight 0.3545, 0.354, 0.3538, 0.3537, 0.3536, 0.3535, and 0.353. The line has theslope corresponding to δ =0.6. dotted straight line has the slope corresponding to δ = 0.6. Inset shows the steady-state density of tumor cells xt as a functionofp. Inthevicinityofthecriticalpoint,xt seemsto decay linearly (β=1). ∆=3. We do notpresentnumericalresults,but simula- tionsshowthatthesameestimationsofcriticalexponents ∆=0 ∆=20 ∆=100 areobtainedfor∆=5and10. Hence,ourresultssuggest that the critical exponents β = 1.0(1) and δ = 0.60(5) are universal most likely for any ∆>0. Let us notice that for increasing∆, the criticalrate pc decreases and the tumor-free phase shrinks. Thus, an- giogenicfactors,thatinourmodelcorrespondtopositive correlationbetweennutrientsupplyandtumorcelldistri- bution(∆>0)makestarvationoftumorstodeathmore FIG.5. Differentmorphologiesofd=2tumorsgrowingfrom difficult. Such a behaviour of our model is very plau- a single tumorcell surrounded byempty lattice sites. Calcu- sible since angiogenesis was ”invented” by tumors just lations weremadeon500x500 lattices andpwas chosen such for this particular reason. For very large ∆, the critical that the steady-state density in all three cases was approxi- rate pc seems to vanish, which means that tumors can mately equal(xt ≈0.05). survive evenunder a very small rate of the nutrient supply. Although the presented model is too simple to de- scribedetailedcomplexityofrealtumors,inouropinion, single tumor cell and all other sites empty. Simulations it qualitatively correctly captures the role of angiogenic were performed for various p and ∆ but such that the factors. steady-statedensity oftumor cells xt wasapproximately Despite possessing a single absorbingstate, our model the same and equal to 0.05. Numerical simulations show mostlikelyforany∆>0doesnotbelongtothedirected- that ∆ strongly affects the shape and dynamics of grow- percolation universality class. The exponent β seems ing tumors (Fig.5). While for small ∆ their shape is to take the classical, mean-field value β = 1 that for highly irregular, for large ∆ nearly circular shape with a d = 1 model is certainly a puzzling result. The a well defined boundary is seen. Apparently, in the for- exponent δ = 0.60(5) is different than the mean-field mer case the growth is very much affected by stochastic value δ = 1 [20] and than the directed-percolation value fluctuations while in the latter their role is diminished. (0.159). Possibly novel critical behaviour of our model Differentmorphologiesofgrowingtumorssuggestthat is confirmed with calculationof other critical exponents, their other characteristics will also strongly depend on butdetailedpresentationofourresultswillbegivenelse- ∆. We measured the tumor survival probability starting where [19, 21]. from the same configurations containing a single tumor We also examined the d = 2 version of our model. cell and monitored whether tumor cells survived until a Obtained results [19] suggest that the values of critical given(large)simulation time. The simulations were per- exponents are very close to those in the d = 1 version. formed for several values of p and ∆ and the survival Thed=2versionofourmodeliscertainlymorerealistic probability Psurv is plotted as a function of the steady- in the context of tumor-growth modeling. We examined state density of tumors xt, that was estimated indepen- their growth starting from a configuration containing a dently with standard steady-state simulations (Fig. 6). 4 1 [2] H. M. Byrne, NatureRev.Cancer 10, 221 (2010). 0.9 [3] E. L. Bearer et al., Cancer Res. 69, 4493 (2009). J. 0.8 T. Oden, A. Hawkins, and S. Prudhomme, Math. Mod. 0.7 Meth. Appl. Sci. 20, 477 (2010). 0.6 [4] J. Paulsonn, Phys. Life Rev. 2, 157 (2005). Psurv 0.5 [5] ESc.i.G2a0b,et1t0a05an(d20E10.)R. agazzini, Math. Mod. Meth. Appl. 0.4 [6] N. Bellomo and M. Delitala, Phys. Life Rev. 5, 183 0.3 (2008). N. Bellomo, et al., Math. Mod. Meth. Appl. Sci. 0.2 ∆=0 20, 1179 (2010). 0.1 ∆∆==130 [7] P. Tracqui, Rep. Prog. Phys.72, 056701 (2009). 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 [8] R. P. Araujo and D. L.McElwain, Bull. Math. Biol. 66, xt 1039 (2004). FIG. 6. Tumor survival probability Psurv as a function of [9] L. Preziosi, Cancer Modelling and Simulation (London: steady-statedensityxt. Initial configuration contained asin- Chapman and Hall, CRC Press, 2003). gletumorcellsurroundedbyemptylatticesites. Calculations [10] M. Alber, M. Kiskowski, J. Glazier, and Y. Jiang, On were made on 200x200 lattice and evolution of tumors was monitored until t = 104. Calculations were made for sev- Cellular Automaton Approaches to Modeling Biological Cells,inIMA134: Mathematicalsystemstheoryinbiol- eralvaluesofpandthecorrespondingsteady-statedensityof ogy, communication, and finance, p.12 (Springer-Verleg, tumors xt was estimated from steady-state simulations, sim- New York 2002). ilarly tothe d=1 case shown in Fig. 1. [11] Y.Jiang,Pjesivac-Grbovic,C.Cantrell,andJ.P.Freyer, Biophys. J. 89, 3884 (2005). Letus notice thatsome generalargumentsdevelopedfor [12] S. Torquato, Phys. Biol. 8 015017 (2011). models with absorbing states suggest that in the active [13] Y. Bar-Yam, Dynamics of complex systems, (Westview Press, 2003). phase, both Psurv and the steady-state order parame- [14] E. D. Schwab and K. J. Pienta, Med. Hypoth. 47, 235 ter should scale with the same critical exponent β [20]. (1996). Thus, one expects an approximately linear dependence [15] A.R.A.AndersonandV.Quaranta,NatureRev.Cancer Psurv ∼ xt for sufficiently small xt. Numerical results 8, 227 (2008). confirm that Psurv increases linearly with xt but they [16] Y. Munsury and T. S. Deisboeck, Modeling tumors as also show a strong dependence on ∆. Indeed, for the complex biosystems: an agent-based approach. In: T. S. same xt, the tumor survival probability for ∆ = 3 and Deisboeck ed. Complex Systems Science in Biomedicine p. 573 (New York,NY:Springer; 2006). 10is muchlargerthanin the ∆=0 case. This is yetan- [17] SincethenumberofsitesN inourmodelwasinsomesim- other indication of the importance of angiogenic effects ulations quite large, N ∼ 106, the standard implemen- on the tumor growth. tations of the roulette-wheel algorithm, based on search In conclusion, we examined a lattice model of tumor techniques, would not be efficient. We used the recently growth where the positive correlation of nutrient supply introduced O(1) implementation of this algorithm (A. with tumor cells distribution mimics angiogenic factors. Lipowski and D. Lipowska, Roulette-wheel selection via Obtained results show that such a correlation shifts the stochastic acceptance, e-print: arXiv:1109.3627). [18] J. Wendykier, A. Lipowski, and A. L. Ferreira, Phys. locationofatumorextinctionandmakestheirstarvation Rev. E. 83, 031904 (2011). to death more difficult. Moreover,it changes the critical [19] A.L.FerreiraandA.Lipowski,inpreparation.Similarly behaviour of the model, even though as a model with a to some related models [18], one can formulate a mean- single absorbing state it should belong to the directed field approximation for this model. The obtained results percolation universality class. Surprisingly, even in the are in good agreement with Monte Carlo, especially in d = 1 version the exponent β seems to take the classi- the d=2 version, and will bereported elsewhere. cal, mean-field value β = 1. However, the exponent δ [20] H.Hinrichsen,Adv.Phys.49,815(2000).G.O´dor,Rev. Mod. Phys.76, 663 (2004). describing the time decayofthe order parameteratcrit- [21] We checked that decay of the order parameter at crit- icalitytakesa non-classicalvalue 0.60(5)andmostlikely icality is described by nearly the same exponent (δ = the model represent a novel critical behaviour. Field- 0.60(5)) thatdoesnot dependon theclassof initialcon- theory methods were applied to various models with ab- figurations. In particular, we made simulations starting sorbingstates[22]andonecanhopethatalsointhiscase from configurations containing various concentrations of they could provide valuable insight. In the d = 2 ver- tumor cells and nutrients. [22] O. Al Hammal, H. Chaté, I. Dornic, and M. A. Munõz, sionsuchacorrelationaffectsthemorphologyofgrowing Phys. Rev. Lett. 94, 230601 (2005). M. Mobilia, I. T. tumors and substantially raises their survival probabil- Georgiev, and U. C. Täuber, J. Stat. Phys. 128, 447 ity. From the statistical mechanics perspective it would (2007). be desirable to examine the behaviour of our model for −1 < ∆ < 0 as well as in the limit ∆ → ∞, where the modelsimplifiesandnutrientsupplyisforbiddentooccur on empty sites. [1] F. Michor, Y. Iwasa, and M. A. Nowak, Nature Rev. Cancer 4, 197 (2004).