Received 14 April 2015; accepted 14 June 2015; published 17 June 2015

ABSTRACT

In an environment that is neither static nor in equilibrium, but is dynamic and changing, the kinetics of the reactions that cause the growth of a tumor, which depend on the state of the evolving environment, cannot be parametrized in terms of constant rates. We propose a simple model for describing the growth on an untreated tumor in such environments, which is characterized by a minimal number of parameters and is generalizable to include the effects of various types of therapies. In the simplest version that we consider here, it consists of a linear equation with a time-dependent growth rate, which we interpret as the coupling of the system with a dynamic environment. A complete solution is given in terms of the integral of the growth rate. The essential features of the general solution are illustrated with a few examples, and comparison is made with the models that have been proposed to describe recent data.

Keywords:

Tumor Growth, Mathematical Modeling, Linear Models, Dynamic Environment, Minimal Parametrizations

1. Introduction

Many applications of the mathematical modeling of tumor growth as a function of time have been based on the Gompertz equation [1] -[3] and a modified version of it [4] , power law equations [4] [5] , and generalizations of the logistic equation [6] [7] . Two common features of these models are that their nonlinear differential equations have solutions that reproduce the S-shaped form of the data for untreated tumors, which represents the inhibition of growth towards an asymptotic limit, and that they contain several free parameters to be adjusted according to experimental data. During the last two decades, an immense amount of effort was devoted to include, for example, the effects on tumor growth due to necrosis, apoptosis, angiogenesis and tumor geometry [8] . Many references to previous work can be traced back from the above papers and from [9] . In addition, including additional differential equations to model therapy have served as a guide to improve the effectiveness of treatment. Therefore, the common strategy along these lines of mathematical modeling of tumors consists in starting with a basic model represented by a nonlinear differential equation, such that by including additional terms in the original equation and introducing other differential equations, the different characteristics of the tumor and effects of therapy can be studied.

However, some of the models consist on generalizing another model by modifying a parameter in a way that seems to be motivated mostly by the fact of being able to fit certain data. For example, the so-called Generalized Logistic model is the result of modifying the logistic equation by inserting a term that contains an arbitrary power of the growth function. In a recent work [7] , the authors find that this model is the best that fits their data, with taken to be a rather arbitrary noninteger number. In addition, in this model and the ones mentioned above, the only solution is the trivial solution if the initial population value is specified as being zero. These nonlinear models may be of course appropriate for particular purposes. However the fact that they are not motivated by more fundamental principles makes it difficult to establish a connection between the observable effects that those models attempt to describe and more fundamental understandings of the mechanism of growth.

On the other hand, a few models of tumor growth that have been recently proposed, are based on assumptions and considerations that take into account basic physical principles, such as energy conservation, scaling, a fluctuating growth rate [10] -[13] , and/or other features specific of the system under consideration. A succinct discussion about cancer modeling that offers the author's insight into the various points of view regarding the mathematical modeling of tumor growth, together with a representative list of references, can be found in [14] .

The main objective of this article is to present a basic model of growth of an untreated tumor in an environment that is not static, nor in equilibrium, but is dynamic and changing. In such a case the kinetics of the reactions that cause the growth of the tumor, which depend on the state of the evolving environment, cannot be parametrized in terms of constant rates. The model borrows some ideas from those mentioned above, contains a minimal number of parameters, and can be generalized to include other effects, such as those of therapies. In its simplest version, which we consider in some detail here, it consists of a linear first-order differential equation with a time-dependent growth rate, which we interpret as being due to the coupling of the system, the tumor in our case, with the dynamic environment. As argued in [14] , simple models can give insights about the workings of detailed and more comprehensive models, and it is in this spirit that we present this class of models here.

As we show, the solution, which can be expressed in terms of an integral of the growth rate, reproduces the usual S-shaped form that vanishes at and eventually reaches an asymptotic form at large t. For illustrative purposes, we consider in Section 3 (Example 3), a specific application to some recent data reported in [7] , and we show that our model can reproduce remarkably well the numerical fit that the authors make of their data. However, in contrast to the mathematical models mentioned above, the model proposed here contains a minimal number of parameters which have a concrete and well-defined meaning and are in principle determined if the interactions that govern the underlying microscopic mechanism of growth are known. We believe that this model and the reasoning leading to it could in turn shed light on such mechanisms, and that it provides a firm footing for pursuing and extending this approach to consider in a simpler way a mathematical modeling of therapy and the case of multicellular populations.

The rest of this paper is organized as follows. In Section 2 we write down the model that we propose, we discuss the framework that motivates it, the assumptions and idealizations involved, and the interpretation of the parameters that appear in it. In Section 3 the general solution to the equation is displayed, and it is illustrated by considering various specific examples that can be of practical use. In particular, here we confront the model with the data presented in [7] to illustrate how well it can fit that data with just one unknown parameter in the growth rate coupling. Finally, Section 4 contains our conclusions. As we discuss in that section in more detail, the approach that we propose here has several attractive features that pave the way for further work. In particular, these linear models contain a minimal number of parameters. Those parameters are related to the interactions that govern the growth mechanisms, and therefore the model can in turn shed light on such mechanisms. In addition, while in this paper we have restricted ourselves to the growth of an untreated tumor, similar ideas can be applied to include the effects of therapy.

2. The Model

Our starting point is the equation that describes the growth of a population that is sustained by an environment which is not necessarily static. In situations in which the population function grows up to a certain asymptotic limit, we assert that the rate of change of the population is proportional to the difference between the actual value of the population and its asymptotic limit, where the proportionality factor, as well as, depend on the dynamical variables that describe the environment. In a changing environment, which includes the case in which the environment may be affected by the population itself, and change with time and the equation is

(1)

In situations in which the environment is static, e.g., large enough that it is not affected by the population, and are constants and this equation is reminiscent of Newton’s cooling law which states that the rate of change of the temperature of a system is proportional to the deviation of the system’s temperature from the temperature of its environment. In our context it may be possible and useful to give a motivation and justification of Equation (1) in terms of more basic principles as follows.

The problem of the time evolution of the population of a given specie appears in many physical contexts. Two basic principles that guide the development of a population are the Master Equation and the Principle of Detailed Balance. The master equation takes the form

(2)

where W depends on f itself and the other variables that describe the rest of the system with which the population can interact. W is decomposed into a series of terms, each of which represents the contribution due to a particular process that causes the population to change. The principle of detailed balance states that there is a precise relation between the so-called direct process and its inverse.

For example, let us consider a process in which only one cell participates and let us denote such process in symbols by

(3)

where stands for a member of the population (a cell) and X stands for a different object. In the direct process, indicated by the right-pointing arrow, a cell disappears into X, while in the inverse process, indicated by the left-pointing arrow, the reverse is true. Then, denoting by and the rates for the direct and inverse processes, respectively, their contribution to W is

(4)

where

(5)

Similar equations also describe the kinetic approach to equilibrium of systems that are put in contact with a reservoir. In such cases, which are governed by physical kinetic equations, the Principle of Detailed Balance can be applied to establish a relationship between and. In our case we can envisage that, in principle, both could be calculated if the interaction between the cells with their surroundings and among themselves were known. However, we have at present no such theory of these interactions. Thus, we leave as an unknown parameter with the property that is a positive quantity.

The procedure outlined above for the case of single cell processes could be generalized to more complicated ones. For example, consider the processes in which two cells participate, which we denote in symbols by

(6)

Because the direct process involves two cells, its rate is proportional to. By the same reasoning that lead us to write Equation (4), the contribution from these processes to W is of the form

(7)

where characterizes the rate for the process to occur.

As a typical rule in those contexts in which these equations have already been applied, the processes in which more than two members participate are rare and not important. Therefore, we are tempted to state that the master equation

(8)

is a good starting point for further exploration of these ideas in the present context as well.

We are interested in the cases in which the environment can change, due to external influences or by its interaction with the population itself. For us this means that the parameters, and that appear in Equation (8), which depend on the state of the environment, change with time. In the absence of a dynamical theory of the interactions of the cells, all we can do is promote them to be functions of time.

In the present paper, we will restrict ourselves to the linear term only, as written in Equation (1). The assumption behind this approximation is that the process in which the cells participate in pairs are rare compared to those in which only one cell participates. Should this linear approximation prove to be inadequate, it could indicate that the pair interactions are important and the quadratic terms in Equation (8) should be taken into account. Overall, we believe that this approach provides a framework for carrying a systematic analysis, based on incremental approximations, on a firm footing and in an organized fashion.

3. Solution

3.1. Static Environment

When and are constants, Equation (1) has the simple solution

(9)

where is the initial population, which can of course be taken to be zero. However, notice that the population reaches the saturation limit independently of the initial value. This contrasts with the solution of the GL model [7] , which requires a non-zero value or otherwise the solution is the trivial solution.

3.2. Dynamic Environment

We consider the case in which the environment can change due to external influences, and therefore and in Equation (1) are functions of time. In this case, Equation (1) is conveniently solved by the Green function method,

(10)

where G satisfies

(11)

with being the Dirac delta function. Furthermore, in Equation (10) is a solution to the homogeneous equation, chosen such that f satisfies the initial condition

(12)

A suitable Green function for Equation (1) is

(13)

where is the unit step function, and is the solution to the homogeneous equation

(14)

satisfying the condition. The function is then uniquely determined as

(15)

and therefore the solution for is given by

(16)

Needless to say, if and are assumed to be constant, then Equation (16) reduces to the solution given in Equation (9). But for any imaginable functions and that could be used to parametrize the changing environment, Equation (16) readily provides the complete solution in terms of two integrals.

A particularly simple form of the solution is obtained when is a constant. In this case, remembering that g satisfies Equation (14), Equation (16) is trivially integrated and yields

(17)

where we have used.

3.3. Examples

In order to illustrate some general features of the solution, we will consider below various specific cases.

3.3.1. Example 1

Let us assume that is a constant, while varies as some (integer) power of t; i.e.,

18)

where a is positive constant and n is a positive integer. First, from Equation (15),

(19)

where we have defined

(20)

for simplicity of the notation. The solution obtained from Equation (17) is then

(21)

In Figure 1 we plot the function, in arbitrary time units (i.e., setting), for the values of the exponent, and taking.

3.3.2. Example 2

As before, we assume that is a constant, but now take as a combination of two monomials; i.e.,

(22)

where are positive integers and are positive constants. Following the same steps as above, the solution obtained from Equation (17) is given by

(23)

where is defined as in Equation (20) and in analogous fashion

(24)

In Figure 2 we plot the function, for various values of n, m and the ratio. We stress that, apart from the initial and saturation values and, respectively, the only unknown and adjustable parameters of these models are the constants that appear in the ansatz for the growth rate function, e.g. the constants a and b defined in Equations (18) and (22).

3.3.3. Example 3

In order to make contact with experimental data, here we confront the model we have presented with the numerical results obtained with the model considered in [7] . We emphasize that we make this comparison to illustrate how well our model can fit experimental data and with no other purpose. Although the model of [7] belongs to the class of models that we have criticized in the Introduction as having some unphysical features, the fact that it fits the data well serves as a benchmark against which to compare numerically our model, or any other model for that matter.

The model of [7] is defined by the equation

(25)

where denotes the growth function, is the effective growth rate constant and K is called the carrying capacity. The solution is written in the form

(26)

where

(27)

Therefore, K determines the saturation limit, while is the initial value. Equation (26) reproduces the Gompertz function for, or the logistic function for, and in this sense it can be considered to be the most general function of this kind that possesses the appropriate asymptotic behavior and shape consistent with the observed data. A common characteristic of these models is that the differential equation, Equation (25), does not allow a zero population value at.

In order to fit their data, the authors of [7] find

(28)

In Figure 3, the solid line represents the function y for the values of the parameters given in Equation (28), and in the same units as K. The broken line is the function f, where we have taken to be given by

(29)

Which has the solution as given in Equation (21), with, namely

(30)

For the plot in Figure 3, we have set and equal to and K respectively, and we have taken

(31)

Obviously, the function will be able to fit the data equally well, with the bonus that, apart from the initial and saturation values and, respectively, the only unknown and adjustable parameter of the solution in our model is the constant that appears in Equation (29). Moreover, in contrast with the adjustable parameters that enter in Equation (26), the parameter characterizes the function which has a definite physical interpretation in terms of the growth rate.

3.4. Generalizations

By inspection, Equation (23) can be generalized in an obvious way to the case in which is a power series

in t. Although the models considered above, with taken to be a constant, already provide a flexible and rich structure as far as their ability to fit the phenomenological data is concerned, it is possible to consider the cases in which is a function of t. In these cases it is not possible to give a closed expression for the solution given in Equation (16), in general. However, in other physical contexts where similar equations and situations arise [15] , some appropriate approximation methods have been used which could be employed in these cases as well.

4. Conclusions and Outlook

We have presented a model for the growth of an untreated tumor in an environment that is dynamic and not in equilibrium. In such environments, the equations that describe the kinetics of the reactions that cause the growth of a tumor, which depend on the state of the evolving environment, cannot be parametrized in terms of constant rates. The model we propose to study these cases is based on ideas borrowed from models used in some physical contexts, together with plausible assumptions and idealizations that pertain to the application at hand here. In the simplest version, the model consists of a linear first-order differential equation with a time-dependent growth rate, which we interpret as being due to the coupling of the system, the tumor in our case, with the dynamic environment. As shown in Section 3, for that case a complete solution can be readily obtained in terms of an integral of the growth rate function. The solutions were given explicitly for a few sample cases, and they exhibit the known characteristic features of tumor growth. Moreover, the model was compared with the mathematical models that have been used to describe recent data, and it was shown that it can fit the data equally well, but without the unphysical features that those models have.

The approach that we have followed is fruitful in several ways. Firstly, the model contains a minimal number of parameters, which have a well-defined meaning, and are in principle determined by the interactions that govern the underlying microscopic mechanism. Secondly, this model could in turn shed light on such mechanisms, thereby providing a basis for pursuing this line of work. Thirdly, the model can be extended beyond the linear approximation that we have used if the terms that we have neglected are believed to be important in a particular system and the quadratic terms in Equation (8) should be taken into account. While in this paper we have restricted ourselves to treat the growth of an untreated tumor, similar ideas can be applied to include the effects of therapy.

References