Chaotic Dynamics of Tumor Growth and Regeneration*

 

Ceferino Obcemea

Memorial Sloan-Kettering Cancer Center

New York, NY

 

1.0 Motivation:

In the analysis of dose-response of tumor tissues to radiation, there have been a flurry of discussions on the need to account for observed tumor regeneration during and/or immediately after dose delivery. This tumor re-growth can at times be dramatic as to alter the intended local tumor control or the patient’s long-term prognosis [1].

Currently, experimental cell-survival curves which form the basis for dose-fractionation schemes could not account for this phenomenon either because the total time duration of observation is too short (typically a few days or weeks) or that transients in the observed data are always averaged and smoothed out for curve fitting.  In addition, the time factor in fractionation schemes like the a/b model is usually suppressed. It is now clear however that during a typical dose fractionation schedule, rapid tumor proliferation could undermine the treatment success: an accelerated clonogen repopulation may lead to post-treatment recurrence or insufficient in-situ tumor ablation.

The various attempts to model tumor growth during the treatment time interval as well as during its unperturbed state have been variations of density-limited kinetic equations. Some of these are the exponential [2], logistic [3], Gompertz [4], Gomp-ex [2], and von Bertalanffy [5] models which have been proposed to model the growth kinetics giving the time behaviour of  intrinsic tumor growth  and its dose-response curve. The justification for these equations depends on how well they fit the survival curves of the particular tumor system under study. The exponential growth is less general as it holds only for short intervals; the Gompertz and the logistic describe more realistic kinetics even for unperturbed tumor systems.

The Gompertz and logistic equations and their variants belong to a wide class of nonlinear differential equations describing density-limited growth. However, we shall see that the locus of these equations could only give rise to a sigmoidal time plot: the tumor population could initially grow very fast, then decelerate after a time lapse and eventually plataeu out to an asymptotic limit. The plot of the growth rate versus population would trace a one-hump graph where the growth rate is explosive when the population is small, reaches a maximum at some intermediate population size and then goes down to eventual zero, as the population increases further, past the critical size. This behaviour seem to realistically model the tumor growth past a critical volume, parts of it are denied nutritive access from the vascular supply lines, giving rise to a necrotic core. Its growth then begins to slow down till it reaches a "quiescent" equilibrium size. Unless further angiogenesis [6 ] creates new supply lines so that it could grow again to the next size-threshold.

The sigmoidal and one-hump functions appear to be nice, simple and well-behaved functions: one or two parameters of the equation could be adjusted to change the steepness of the sigmoid plot or the amplitude or the skewness of the hump. They are thus ideal for curve-fitting, least-squares optimization or linear regression analysis vis-à-vis an experimental tumor growth curve or when coupled to the a/b model, against a dose-survival curve.

This simplicity could however be disarming and may prove illusory. The purpose of this note is to bring notice to the fact that due to its nonlinearity, this class of functions actually exhibit very complicated behaviour. To wit, as the value of their parameter e.g. proliferation rate coefficient increases, the approach to the asymptotic limit instead of being smooth becomes wildly jagged, then the asymptote itself bifurcates and is no longer single-valued and the population simply oscillates between two population points. At a higher value of the parameter, the population loses any concept of a fixed asymptote; the population plunges into limit cycles, going back into some set population points after four, eight and higher periods. Finally, at a critical value of the parameter, the oscillations become unyielding, the possible values of the population become manifold; a time-behaviour that has been termed "chaotic".

One may argue that this chaotic behaviour is far and between and could only happen on rather rare exotic occasions. And yet such phenomena have already been demonstrated to occur in biological systems such as circadian rhythms, EEG wave forms, ion transport across cellular gap junctions, chaotic oscillations in tumor model systems [7] as well as population dynamics of many predator/prey species [8].

Where is the origin of the chaotic dynamics of tumor growth? Why is chaos not manifested in the time-integrated form of the logistic or Gompertz equations? The biological basis of this dynamics comes from the fact that tumor cells grow as with normal ones, in discrete time intervals. Growth happens with defined cell cycles with characteristic doubling times from the parent cell to the progeny after mitosis. Hence, the underlying equations should be discrete maps and not continuous differential equations. In the continuous case, time integration between two end-points smooths out oscillatory behaviour, which may be present. Discretizing this integration interval would recover these oscillations. In fact, the fortuitous discovery of chaotic behaviour in ecological models [8] came from the numerical iterates of the discrete version of the underlying differential equation. We thus convert the differential equation over an interval into a difference equation. Rather than a one-time integration between end-points of the interval, one could follow the behaviour of the solution in discrete steps, going from one iterate to another until we exhaust the interval.

What is the importance of knowing this complicated growth dynamics? For one, it gives the caveat against adjusting the parameters of the growth function simply to curve-fit the dose–survival data. Large number of parameters may give a good fit but may as well give unstable solutions as the order parameter increases: the instability being unnoticed because most observation times e.g. in experimental tumors in-vitro or in animal models are too short. Also, the error bars in many experimental growth curves may prove not to be "errors" or "uncertainty" but real data point dispersions as some tumor model systems exhibit different growth patterns when experiments get repeated with small variations of physiological conditions.

More importantly, nonlinear dynamics may help us understand the problem of accelerated tumor clonogen repopulation [9]. It may be that the change in the doubling times between pre- and post-treatment of target tissues could catapult the clonogen population to the chaotic regime or at least to limit cycles of higher periods. This change in the doubling times could be drastic, being 14-fold for example in head and neck carcinomas.

 

2.0 Continuous-time growth models and its limitation

The logistic and Gompertz equations are two popular models currently being employed to describe tumor growth kinetics. They come from models of population dynamics, which reflect inhibitory effect of population density:

dN/dt = aN –bN2 (Logistic)                                                                       (2.1)

dN/dt = aN – bN lnN (Gompertz)                                                                                            (2.2)

 

where a and b are the parameters of the equation.

Thus, instead of a runaway Malthusian growth, either the square or the logarithm of the population forces the growth to plataeu out. We focus on the logistic equation but the result hold as well to the general class of one-hump functions. In the logistic equation, a and b are the proliferation rate and density coefficients: a gives an exponential increase when the population is small and b damps out the growth rate as the second term of (1.1) predominates when N is large. On integration,

Nt = aN0 / [bN0 + (a –n N0) exp (-at) ]                                                                    (2.3)

so that as t increases, Nt asymptotically goes to the limiting value: a/b., often called the carrying capacity of the population. Graphically, the time plot of Nt is a sigmoid: as Nt becomes large, the logistic support each member of the population receives becomes more scarce and the growth rate slows down. What is the maximum growth rate that the population could achieve?

d/dN [ (dN/dt)] = a – 2bN = 0                                                                                  (2.4)

(dN/dt) at Nmax = a2/4b                                                                                           (2.5)

(Fishermen have already known this result since Volterra’s analysis: what maximum number of fish could one catch per unit time while ensuring that the fish population still remain viable?). In the context of tumor population, this tells us that any cellular assault due to regimens of radiation or cytotoxic drugs that can put a decrement on the tumor population equal to (a2/4b) would control the growth. At a rate greater than (a2/4b), this would, in time, render the tumor population extinct. Plotting the growth rate with respect to the population level N gives a one-hump function with the value of (a2/4b) at maximum height. Different values of a and b give a family of curves of varying steepness of the sigmoid and amplitude of the hump.

Note that no other time behaviour is manifest in this curve. Observationally, indeed tumor size can approach an asymptotic value, growing to a certain volume where it seems to max out, i.e. becomes quiescent. However, it is also known that tumor often spontaneously regress [10], especially when they are very small, possibly unable to outwit immune surveillance. Hence, tumors also become extinct as a matter of time-course. In addition, it is also known that it can grow, regress, re-grow in an oscillatory manner. And then finally, growth can sometimes remain rampant without obvious limit, until the host itself perishes from the sheer tumor burden. These time-behaviours are never accounted for by current continuous-time tumor growth models.

 

3.0 Discrete Maps as realistic tumor growth model

As we stated before, tumor growth kinetics should be viewed in discrete times, since growth occurs in distinct cell cycle with characteristic doubling times. Hence, the more realistic version of (2.1) should be its discrete analogue. Even in the context of dose-response of tumor tissues, dose is always delivered in discrete amounts and/or discrete time intervals. The nonlinear difference analogue of (2.1) is then

 

Nt+1 = Nt + a Nt – b Nt 2                                                                               (3.1)

where Nt+1 is the population after one generation. It is clear that Nt = a/b is the asymptote i.e. the next generation Nt+1  will never surpass the previous one when Nt reaches the value of a/b. If we change variables to U= N(b/a), α = a +1, then (2.1) becomes the familiar logistic map:

Ut+1 = α Ut (1- Ut  )                                                                                      (3.2)

It well known that the analysis of the iterates of this map gives rise the various time-behaviours described above: namely spontaneous extinction or regression, oscillatory growth into various limit cycles and finally a period-doubling cascade to chaos.

What is also interesting is the clinical origin of the parameters a and b . The advances in understanding the genetic machinery of tumor growth and molecular mechanism in which various cytokines and growth factors could trigger this machinery give indication on how these parameters change. a is an intrinsic growth rate that can be modified by defects in cellular genetic machinery e.g. p53 gene deficient or mutated [11], by stimulatory growth factors such as TGFs, EGFs, angiogenic growth factors and by various cytokines [12]. The carrying capacity of the tumor system for example can be dramatically increased by the angiogenic or p53 mechanism. b is an inhibitory parameter that can be modified by inhibitory growth factors such as TGFb , NGFs and /or cytokine that stimulates the immune response such as the interleukins and interferons.

The analysis of the bifurcation structure of the discrete maps such as above, together with the clinical elucidation of the system parameters a and b provides a promising approach towards understanding this difficult dynamics of tumor growth.

 

References

  1. J. Cox,T Pajak, V. Marcial, L. Coia, M. Mohhiuddin, K. Fu, H. Selim, R. Byhardt, P. Rubin, H. Ortiz and L. Martin (1992), Cancer, 69, 2744.
  2. G. G. Steele, (1977), Growth Kinetics of Tumors, Clarendon Press (Oxford).
  3. A. Laird, (1964), Br. J. Cancer, 18, 490.
  4. J. Speer, V. Petrosky, M. Retsky and R. Wardwell, (1984), Cancer Research, 44, 4124.
  5. Z. Bajzer, M Marucic, S. Vuk-Pavlovic, (1996), Math. Comp. Modelling, 23, 31.
  6. J. Folkman, (1995), J. Molecular Med., 1, 120.
  7. E. Posadas, S. Criley, and D. Coffey, (1996) Cancer Research, 56, 3682.
  8. R. May, (1974), Science, 186, 645.
  9. H. Withers, J. Taylor and B. Maciejewski, (1988), Acta Oncologica, 27, 131.
  10. G. Challis and H. Stam, 1990, Acta Oncologica, 29, 545.
  11. T. Graeber, C. Osmanian, D. Housman, C. Koch, S. Lowe and A. Giaccia, (1996),Nature, 379, 88.
  12. M. Sporn and A. Roberts, (1985), Nature 313, 745.

 

*Paper given at the 3rd Int'l. Conference on Chaos and Complex Systems, Nashua, NH May 21-26, 2000.

  To appear in Proceedings of the Third ICCS, Perseus Books, Boston 2001.

  Manuscript published in the online journal: Interjournal on chaos and complexity: www.interjournal.org,

  (paper number: 424).