Linear model of tumor growth in a changing environment

We propose a model for describing the growth on an untreated tumor, which is characterized in a simple way by a minimal number of parameters with a well-defined physical interpretation. The model is motivated by invoking the Master Equation and the Principle of Detailed Balance in the present context, and it is easily generalizable to include the effects of various types of therapies. In the simplest version that we consider here, it leads to a linear equation that describes the population growth in a dynamic environment, for which a complete solution can be given in terms of the integral of the growth rate. The essential features of the general solution for this case are illustrated with a few examples.


Introduction
In the last few decades, the mathematical modeling of tumor growth as a function of time has been based mainly in applications of the Gompertz equation [1,2,3] and a modified version of it [4], power law equations [4,5], and generalizations of the logistic equation [6,7]. These mathematical models helped to provide an understanding of tumor growth as a more organized phenomenon than previously * nieves@ltp.upr.clu.edu † ubriaco@ltp.upr.clu.edu thought. In particular, solutions of the Gompertz and the modified logistic equations exhibit an Sshaped form which represents inhibition of growth to an asymptotic limit. In addition, their applications to experimental data, with the help of additional differential equations to model therapy, has served as a guide to improve the effectiveness of treatment. However, a preference for one model or the other has been based exclusively on their adequacy to fit experimental data. Similarly, some of the models consist on generalizing previous ones by modifying a parameter in a way that is not motivated by any fundamental principle, but again to fit some data. For example, the so-called Generalized Logistic model is the result of merely modifying the logistic equation by inserting an arbitrary power which, in order to fit the data, becomes a noninteger number.
While this may be appropriate for particular purposes, and for applications to systems that are not governed by the laws of Nature, the fact that the models are neither motivated nor based on fundamental physical principles prevents us from ascribing a physical meaning to the parameters that appear in such models. This makes it more difficult, if not impossible, to establish a connection between the observable effects that can be described on the basis of phenomenological models, with more fundamental explanations and understandings of the mechanism of growth in biological systems, which must be based on detailed microscopic dynamics.
The model is motivated by invoking the Master Equation and the Principle of Detailed Balance in the present context. In the simplest version that we consider here, it leads to a linear equation that describes the population growth in a dynamic environment, The main objective of this article is to present a model that is inspired by fundamental physical principles, such as the Master Equation and the Principle of Detailed Balance, and is partly motivated by analogous equations that apply to a variety of physical systems. The model is based on a linear differential equation with a time dependent coupling. It shares the same asymptotic behavior of the GL model, giving an S-shaped function for tumor growth that vanishes at t = 0 and eventually reaches an asymptotic value at large t. However, the model proposed here contains a minimal number of parameters, which have a concrete and well-defined meaning and are in principle determined and calculable if the interactions that govern the underlying microscopic mechanism of growth are known. By the same token, the use of this model should in turn shed light on such mechanisms, thereby providing a firm footing for pursuing and extending such approaches.
The rest of this paper is organized as follows. In Section 2, for the purpose of comparing our proposed differential equation with the previously proposed models, we briefly recall the basic aspects and shortcomings of them. In Section 3 we discuss the theoretical framework that motivates the model that we propose, the assumptions and idealizations involved, and the interpretation of the parameters that appear in it. In Section 4 the general solution to the equation is displayed, and it is illustrated by considering various specific examples that can be of practical use. Finally, Section 5 contains our conclusions.

Previous models
We will denote by f (t) the growth function, and for simplicity of the notation omit the argument t except when necessary to avoid confusion. We envisage f to be the number of cells at a particular instant of time.
The Gompertz equation can be written in the form where α and β are two parameters. The solution is given by with where f 0 is the initial value of f . The asymptotic value is determined as A particular feature of Eq. (1) is that it is not defined for f = 0, and therefore it does not allow the initial value f 0 = 0, as Eq. (3) reveals, thus ruling out its application to data with rather small initial population value. A way to overcome this difficulty is to modify Eq. (1) by introducing an arbitrary parameter ǫ > 0 such that the new differential equation reads [4] which has a non trivial solution given by where While this procedure overcomes the above-mentioned difficulty with the Gompertz solution, a physical motivation or justification is lacking. Another approach [5] has been to consider the power law differential equation For β = 1 it has the solution which leads to linear or exponential growth for α = 0 or α = 1, respectively. Certainly, power law growth is unconstrained, and its behavior is radically different from that given by Eq. (1).
Another way to obtain asymptotic behavior in tumor growth is to modify the logistic equation by introducing an arbitrary power ǫ > 0, leading to the differential equation whose solution is given by where Eq. (11) reproduces the Gompertz function for ǫ = 0, but it also has the limitation that the differential equation Eq. (10) does not allow a zero value at t = 0. Therefore, the solution of Eq. (10) can be considered to be the most general function that possesses asymptotic behavior and a shape consistent with experimental data. However, the fact that Eq. (10) can be neither motivated nor understood on the basis of some of fundamental physical principles is in our opinion an important limitation to further understanding the mechanisms that drive the growth in these biological systems, on the basis of such equations. The main reason for this lack of insight is due to the fact that Eq. (10) is the result of merely modifying the logistic equation by inserting an arbitrary power, that in order to fit the data becomes a noninteger number. Therefore, without any other justification it is not possible in this and the previously discussed models to ascribe a deeper physical significance to the parameters of the model.

The Model
Our starting point is the equation that describes the growth of a population that is sustained by an environment. We assume that in such situations the population grows up to a certain saturation limit f s , and that the environment is large enough such that it is not affected by the population itself. Under such conditions, we assert that the rate of change of the population is proportional to the difference between actual value of the population and its saturation limit.
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 is possible and useful to give a motivation and justification in terms of more basic principles as follows.

The Master Equation
The problem of the time evolution of the population of a given specie appears in many physical contexts. For example, in the astrophysical context of the Early Universe, one analogous problem is the determination of the abundance of the various atomic elements and how they form [8].
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 φ ↔ X , 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 γ d and γ i the rates for the direct and inverse processes, respectively, their contribution to W is written in the form where 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 implies a fundamental relation where ∆E and T are identified with the change in energy of the system and the temperature of the environment, respectively. In our case, in principle both γ d,i 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 can be generalized to more complicated ones. For example, consider the processes in which two cells participate, which we denote in symbols by Because the direct process involves two cells, its rate is proportional to f 2 . By the same reasoning that lead us to write Eq. (16), the contribution from these processes to W is of the form 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 is a good starting point for further exploration of these ideas in the present context as well.
In the present paper, we will restrict ourselves to the linear term only, as written in Eq. (13). 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 Eq. (21) should be taken into account. Overall, this approach provides a framework for carrying a systematic analysis, based on incremental approximations, on a firm footing and in an organized fashion.

Static environment
When γ is a constant, Eq. (13) has the simple solution where f 0 is the initial population, which can of course be taken to be zero. However, notice that the population reaches the saturation limit f s independently of the initial value f 0 . This contrasts with the solution of the GL model [7], which requires a non-zero value f 0 or otherwise the solution is the trivial solution f (t) = 0.

Dynamic environment
We consider the case in which the environment can change due to external influences. For us this means that the parameters γ and f s that appear in the model equation, both of 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 γ and f s to be functions of time. Therefore, our basic equation becomes This equation is conveniently solved by the Green function method, where G satisfies f h is a solution to the homogeneous equation such that f satisfies the initial condition f (0) = f 0 . In Eq. (25), δ(x) stands for the Dirac delta function. A suitable Green function for Eq. (23) is where g(t, t ′ ) is the solution to the homogeneous equation satisfying the condition g(t ′ , t ′ ) = 1, and θ(x) is the unit step function. The function g(t, t ′ ) is then uniquely determined as and therefore the solution for f (t) is given by Needless to say, if γ and f s are assumed to be constant, then Eq. (28) reduces to the solution given in Eq. (22). But for any imaginable functions γ(t) and f s (t) that could be used to parametrize the changing environment, Eq. (28) readily provides the complete solution in terms of two integrals. A particularly simple form of the solution is obtained in the case in which f s is a constant. In this case, noting from Eq. (27) that then Eq. (28) yields where we have used g(t, t) = 1.

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

Example 1
Let us assume that f s is a constant, while γ varies as some (integer) power of t; i.e., where a is positive constant and n is a positive integer. First, from Eq. (27), where we have defined for simplicity of the notation. The solution obtained from Eq. (30) is then In Fig. 1 we plot the function f /f s , in arbitrary time units (i.e., setting γ 0 = 1), for the values of the exponent n = 1, 2, 3, 4, 5, and taking f 0 = 0.

Example 2
As before, we assume that f s is a constant, but now take γ as a combination of two monomials; i.e., where n, m are positive integers and a, b are positive constants. Following the same steps as above, the obtained from Eq. (30) is given by where γ 0 is defined as in Eq. (33) and in analogous fashion In Fig. 2 we plot the function f /f s , for various values of n, m and the ratio r = γ 0 /γ ′ 0 . We stress that, apart from the initial and saturation values f 0 and f s , respectively, the only unknown and ajustable parameters of these models are the constants that appear in the ansatz for the growth rate function γ(t), e.g. the constants a and b defined in Eqs. (31) and (35).

Generalizations
By inspection, Eq. (36) can be generalized in an obvious way to the case in which γ(t) is a power series in t. Although the models considered above, with f s taken to be a constant, already provide a flexible and rich structure as far as their ability to fit the phenomenlogical data is concerned, it is possible to consider the cases in which f s is a function of t. In these cases it is not possible to give a closed expression for the solution given in Eq. (28), in general. However, in other analogous physical problems where similar situations arise [8], very effective approximation methods have been used which could be employed in these cases as well.

Conclusions and Outlook
In this article we have presented a model for tumor growth, which is based on physical principles on one hand, together with plausible physical assumptions and idealizations on the other. As shown in Section 4, in the simplest version of the model, in which the pairwise interactions between the cells are neglected, the growth equation is linear and a complete solution can be readily given. Moreover, the solutions were explicitly given for a few sample cases, which exhibit the known characteristic features of tumor growth. Thus the approach that we have followed is fruitful in several ways. Firstly, the model contains a minimal number of parameters, which have a concrete and well-defined meaning, and are in principle determined and calculable if the interactions that govern the underlying microscopic mechanism of growth are known. Secondly, by the same token, the use of this model should in turn shed light on such mechanisms, thereby providing a firm basis for pursuing this line of work. Thirdly, the model can be extended beyond the linear approximation that we have used if the pairwise interactions are believed to be important in a particular system and the quadratic terms in Eq. (21) should be taken into account. Our approach provides a framework for taking into account such higher order terms in a systematic fashion. Lastly, while in this paper we have restricted ourselves to treat the growth of an untreated tumor, our work paves the way for applying similar principles and ideas to include the effects of therapy. Work along these lines is in progress.