Abstract

A well-known mathematical model of radially symmetric tumour growth is revisited in the present work. Under this aim, a cancerous spherical mass lying in a finite concentric nutritive surrounding is considered. The host spherical shell provides the tumor with vital nutrients, receives the debris of the necrotic cancer cells, and also transmits to the tumour the pressure imposed on its exterior boundary. We focus on studying the type of inhomogeneity that the nutrient supply and the pressure field imposed on the host exterior boundary, can exhibit in order for the spherical structure to be supported. It turns out that, if the imposed fields depart from being homogeneous, only a special type of interrelated inhomogeneity between nutrient and pressure can secure the spherical growth. The work includes an analytic derivation of the related boundary value problems based on physical conservation laws and their analytical treatment. Implementations in cases of special physical interest are examined, and also existing homogeneous results from the literature are fully recovered.

1. Introduction

Mathematical modelling of cancer tumour growth helps in understanding the mechanisms underlying the phenomenon in many ways. The main idea is that modelling physical hypotheses in mathematical terms, a biological phenomenon, such as tumour growth, is approached by a mathematical problem. Analytic and numerical or hybrid methods, applied to the problem studied, lead to conclusions on the solution of the problem. The conclusions are interpreted in biological terms that are subject to evaluation with respect to experimental evidence. Thus, mathematical modelling may potentially offer new perspective in the research directions of biological procedures.

As cancer research and the related technology improved, an increasing amount of experimental data was produced, which made it hard to interpret. At the same time, mathematical models begun to develop at the aid of understanding the crucial parameters involved in tumour’s evolution [1, 2]. In the early years, mathematical models focused mainly on the avascular phase of the tumour’s development. This phase corresponds to the stages right after tumorigenesis and ends with a steady state, where the tumour’s volume gain, due to the new cancerous cells birth, balance the volume loss from the cells’ death and disintegration. In order for the tumour’s evolution to proceed, new phenomena, as angiogenesis, have to take place and the vascular phase may begin. In this phase, the tumour has developed a vascular net around it that provides the tumour cells with limitless nutrient supply and also it permits metastasis. As the present work focuses on the avascular phase on tumour evolution, we avoid presenting details on the proceeding phases, which can be found in many references, see, for example, [1, 2]. Since these next phases include the phenomena, appeared in the avascular phase and other more complicated ones, in depth studying of the avascular phase is well justified.

Both deterministic, numerical and hybrid models have been developed towards the avascular growth study [2]. In the present work, we study a deterministic mathematical model that considers the tumour colony as a continuous medium, which evolves according to mass conservation law, Fick’s diffusion law, and fluid mechanics principles [3]. Deterministic models in tumour growth research have been studied since 1966 when Burton [4] modelled the tumour’s evolution with a Gompertzian curve. Though, the basis for mathematical modelling of avascular tumour growth has been established by the Greenspan works [5, 6] in the 1970s. Greenspan formulated the main procedures of the phenomenon, following existing experimental results, solving boundary value problems with reaction-diffusion partial differential equations in spherical domains. The main parameters that he investigated were the nutrient concentration in the tumour and its surrounding, as well as the concentration of a growth inhibitor factor. In addition, he introduced the idea of a surface tension that secures the solid structure of the tumour. Since then, his ideas were used widely and produced sophisticated models that investigate different aspects of avascular tumour growth. Indicatively, we refer to Adam’s papers [7, 8], to the series of works by Byrne and Chaplain [912], by Ward and King [13, 14], as well as by Mc Elwain and his collaborators [15, 16]. One can have a clear view on relevant references and analytic presentation of each contribution in the extensive reviews [1, 2, 11].

In most works, a spherical tumour structure has been used and different forms of the physical parameters involved have been investigated. Heterogeneity with respect to different cell types has been taken into account [11], and the effect of inhomogeneity of the consumption rates and proliferation rates within the tumour colony has been investigated [15]. However, most works consider homogeneous nutrient supply from an infinite environment whereas this assumption ignores the realistic fact that tumours grow in finite organs with nonhomogeneous oxygen or glucose concentrations in the tumour’s microenvironment. Moreover, there is experimental evidence [17] that a geometrical confinement has significant impact on tumour growth, and to our knowledge there are only few models that take this fact into account and stress out the importance of the tumour-host boundary interactions and the environmental heterogeneity [18, 19]. In the present work, we investigate the way that some key growth parameters interact, so as to permit a spherically symmetric tumour growth. The parameters that we consider are the finite growth environment, the nonhomogeneous pressure that it might induce, and the nonhomogeneous nutrient supply. Although a work on a spherical model is quite trivial after so much previous work on that geometric configuration, the point of view in our paper is to have a close look upon the detailed conditions that allow for such a configuration to occur and so to get the perspective of the conditions that might allow for other, less symmetric, geometries to appear, like the ellipsoidal one. Therefore, the present work should be viewed as a theoretical analysis of the basic mathematical model that provides proof that under the common assumptions the symmetric growth is inevitable and that radial growth is also possible under certain diversions from such assumptions. Moreover, concerning the cells’ motility that results to the tumour’s spatial expansion, we assume the rather logical approach that, additionally to the nutrient gradient and to the pressure gradient [6, 19], we further attribute the cells’ movement to the inhibitor gradient, a parameter that shows up to have a qualitative impact to the final results.

In particular, in this work, we consider a spherical fully developed avascular tumour, meaning that it consists of a necrotic core, a quiescent layer, and a proliferating layer, which are under dynamical growth until the whole structure reaches a steady state. Then, the tumour may stay in a dormant state or it may further develop only by entering in the vascular phase, a case that is not examined in the present study. We further assume that the tumour grows in a finite concentric spherical host environment, modelling the healthy part of the organ in the interior of which the tumour develops. Both tissues, cancerous and healthy, are assumed to be incompressible fluids of different phase. Here, we need to clarify that the tumour’s avascular phase ends at a maximum tumour size, which is significantly smaller than a typical human host organ. Nevertheless, our study aims to reveal rather qualitative aspects than quantitative effects of the tumour growth. One such aspect is concerned with the confinement effects imposed by the growth environment. The effect of the stress, imposed by the surrounding, has been subject of investigation in several works [1, 8, 17, 18]. For example, the effects of a spatially restricted growth have been studied experimentally in [17], where a cancer cell culture was grown inside a very thin cylindrical tube. It turned out that the tube altered the shape and growth dynamics of the cell colony. Similar result was the outcome of the computational work [19]. Therefore, an analytical study that further investigates such dynamic growth inside a finite host environment is justified.

We assume that the tumour receives the nutrient (in general oxygen or glycose) from the surrounding spherical layer with inhomogeneous concentration and also it is affected by a tissue oriented nonhomogeneous pressure field. The proof that either no kind of inhomogeneity or a special kind of inhomogeneity that relates the externally supplied nutrient field with the externally imposed pressure field can be supported by such model in order for its concentric spherical structure to be secured is one of the outcomes of the present work. In the first case, it turns out that a spherical tumour can be developed only under homogeneous nutrient supply and homogeneous pressure field imposed, which confirms all previous works in the field. Nevertheless, the model reveals another option that allows nonhomogeneity in a concentric spherical growth, provided that a special relation holds between the nutrient and the pressure fields, with respect to the cancer cells’ motility parameters. In any case, the concentration profile of both the nutrient and the inhibitor turns out to be strongly affected by the host tissue boundary. Similarly, the pressure distribution throughout the cancer tissue and the host tissue is affected by the host tissue boundary, but the evolution of the tumour boundary is independent of it.

Before we proceed, we present briefly the basic assumptions we make within this work following the framework of [5, 6, 9, 12].

We consider a spherical shell with radius ????, occupied by incompressible noncancerous tissue, which covers a growing spherical tumour colony. We assume that the tumour is avascular and fully developed, while it consists of three concentric spherical domains. The exterior spherical shell O?? is occupied by proliferating cells, the subsequent shell O?? is occupied by quiescent cells, and the interior necrotic sphere O?? includes the dead cells and debris. The crucial physical parameters that govern the tumour growth at every fixed point ??(??) are a nutritive substance supplied externally (usually considered as oxygen or glucose) with concentration ??(??), a growth inhibitor factor with concentration ??(??), which is produced by the necrotic cell debris and the pressure field ??(??) generated by the incompressible finite environment and the growing mass in its interior. The nutrient is supplied into the host tissue with constant concentration ??8(????) and diffuses inwardly with diffusion constant ????. Tumour cells consume nutrient and proliferate, as long as the nutrient concentration is maintained above the critical value ??*??, which is fulfilled in the exterior shell O??. In the inner shell O??, where ??(??) lies in the interval 0<??*??<??<??*??, the tumour cells cannot proliferate and live in a quiescent state. Finally, in O??, the concentration falls below ??*?? and the tumour cells die out of starvation. Moreover, O?? also includes the dead cells due to apoptosis or programmed cell death, which are forwarded there by pressure gradients. The disintegration of the dead cells into simpler compounds produces a growth inhibitor substance in O??, which diffuses outwardly with diffusion constant ????. This inhibitor factor prohibits cell proliferation, unless its concentration falls below a critical value ??*??, which is assumed to hold in O??. As both the nutrient and the inhibitor diffuse inward to or outward from the tumour, respectively, their distribution is described with parabolic partial differential equations, as we will show in the sequel. Though, it is a common assumption [1, 3, 9] that the diffusion time scaling is significantly shorter than the time scaling one is really interested in, which is the tumour volume doubling time scale. Hence, as we demonstrate in Section 2.1.1, the elimination of the time derivative is justified concerning the nutrient and the inhibitor concentration, which are then considered to be in diffusive equilibrium state. The pressure field is generated from the pressure ??8(????), which is imposed by the host boundary and affects the interior pressure following the Young-Laplace law [5, 6] for two face incompressible fluids, together with the pressure imposed by cell proliferation and loss of necrotic mass, due to disintegration. The produced pressure differentials cause passive cell motion opposite to the direction of the pressure gradient. Additionally to the widely used Darcy’s law, which connects the velocity ??(??) of the tumour cells at the exterior tumour surface with the pressure gradient [3, 19], we moreover assume that exterior tumour cells actively move towards the nutrient gradient [18, 20] and also we suggest that they move opposite to the growth inhibitor gradient. This assumption leads to the evolution equation of the tumour’s exterior boundary, which is a highly nonlinear ordinary differential equation and its solution is the goal of this work, the exterior tumour boundary as a function of time. Finally, we assume that, as the tumour grows, it maintains its concentric spherical structure, so all the interfaces of the model are time-dependent functions, connected between each other and the whole model is formulated as a free boundary problem.

In Section 2, we derive the corresponding partial differential equations that describe the nutrient concentration and the inhibitor concentration field, and the pressure field in each region, applying appropriate physical laws. Also, the corresponding boundary conditions are derived in each interface. Finally, the corresponding boundary value problems are stated in the spherical coordinate system. The solution of these boundary value problems is included in Section 3, where the restrictions on the kind of the inhomogeneity that can be supported with the spherical structure are derived. Also, the connection between the inner boundaries and the evolution equation of the exterior tumour boundary is provided. In Section 4, we demonstrate our results by presenting particular implementations regarding the exterior supply of the nutrient and the pressure imposed to the exterior boundary of the host surrounding of the tumour. Also, reduction of the results obtained in the present work to already existing corresponding expressions included in the original works [5, 6, 11] is provided in Section 4. Finally, an outline of our work follows in Section 5.

2. Statement of the Problem

Let (??,??,??) be the spherical coordinates of the point ???R3. Then, each distinct region in the tumour is defined as follows. The necrotic core is defined as O??=?(??,??,??)?R3:0=??<?????,,0=??=??,0=??<2??(2.1) while the quiescent layer, the proliferating layer, and the exterior layer occupied by the noncancerous host tissue are defined as O??=?(??,??,??)?R3:????<??<?????,O,0=??=??,0=??<2????=?(??,??,??)?R3:????<??<?????,O,0=??=??,0=??<2????=?(??,??,??)?R3:????<??<?????,,0=??=??,0=??<2??(2.2) respectively. Each spherical region is bounded externally by the spherical surface ????:??=???? with ??=??,??,??,??.

Figure 1 shows indicatively, with no reference to their real relative size, the aforementioned spherical regions under consideration along with their corresponding boundaries. In reality, the necrotic core is the dominant region of the tumour, while the exterior host layer is much larger than the included tumour regions.

2.1. The Nutrient and the Inhibitor Concentration Problem
2.1.1. Derivation of the Partial Differential Equations

According to the mass conservation law in each of the regions O??, ??=??,??,??,??, the rate of change of the nutrient concentration in each region equals to the normal flow of nutrients ??in(??) inwards to O??, minus the flow outwards, ??out(??), and minus the amount of nutrient that is consumed in the interior of O??, that is, ?O????????(??)?????????=??O??(-??)·??in?(??)????-??O????·??out?(??)????-O??G??????,(2.3) where ??O?? is the boundary of the region O??, ?? is the outward unit normal vector to ??O??, and G?? is the nutrient consumption rate, considered to be constant in O??. Since the nutrient flow is directed monotonically inwards in all spherical layers, the flow inwards takes place in the exterior boundary surface, denoted by ??+??, and the flow outwards is realized from the inner boundary ??-??. Considering also that for the spherical surface the normal direction is the radial one, (2.3) becomes ?O????????(??)?????????=??+??(-??)·??in?(??)????-??-??(-??)·??out?(??)????-O??G??????with??=??,??,??,??.(2.4)

Assuming that the nutrient flow is only diffusive, it follows Fick’s law, which means that it is directed to regions of lower nutrient concentration ??(??)=-???????(??). Then, relationship (2.4) yields ?O????????(??)?????????=??+???????·?????????(??)????-??-???????·?????????(??)????-O??G??????.(2.5)

Applying the divergence theorem in O?? in the right-hand side of (2.5), (2.5) becomes ?O??(((??????(??))/????)-?????????(??)+G??)????=0, which results in ??????(??)????-?????????(??)+G??=0for???O??with??=??,??,??,??.(2.6)

In particular, considering the physical assumptions made in the Introduction, the consumption rate in each region is defined as G??=G??=0, in O?? and O??, respectively, where no consumption of nutrients takes place. In O?? and O??, the consumption rate is assumed to depend on the vital state of the nutrients, which depends on the amount of the available nutrient, while it is defined as G??=G(??*??/??8)and G??=G(??*??/??8), respectively, where G is a proportional constant depending on the species of the host tissue cells, ??*??, ??*?? are the nutrient critical values mentioned in the introduction, and ??8 is a reference value corresponding to the nutrient concentration in the vicinity of the tumour’s host organ.

We turn now to the inhibitor concentration. From the mathematical point of view, the inhibitor concentration problem is very similar to the nutrient concentration problem. The difference lies on the interpretation of the parameters and in their particular expressions. As mentioned in the introduction, the inhibitor is produced inside the necrotic core with production rate ???? and is consumed in the quiescent and proliferating layer with consumption rate ???? with ??=??,??, respectively. Therefore, the inhibitor flow is directed outwards and the corresponding partial differential equations read as follows ??????(??)????-?????????(??)-????=0for???O??,??????(??)????-?????????(??)+????=0for???O??with??=??,??,(2.7) and considering that the inhibitor is neither consumed nor produced in the host tissue, we obtain ??????(??)????-?????????(??)=0for???R3-?O???O???O???.(2.8)

Before we proceed, we will examine the significance of the terms with the time derivatives in the partial differential (2.6)–(2.8). Assuming the reference variables ????, ???? and ??8, we define the dimensionless variables ??=(??/????),??=(??/????)and ????(??)=????(??)/??8for ??=??,??,??,??. Hence, (2.6) becomes ????????????????-???????????+??2????8????G??=0with??=??,??,??,??,(2.9) where ???=??2???, while the constant ??=??2??/(????????) becomes much smaller than one, when the scaling ????=10-2?cm, ????=1day=86400sec, and ????=10-6cm2/sec is used, which are a typical length scale, a typical time scale and a typical diffusion constant for the growth of an avascular tumour [9]. Then, ??=??2??/(????????)=(10-4/(86400·10-6))=1,12·10-3, and the time derivative terms can be neglected. Under the same scaling, (2.7)-(2.8) can also be considered as elliptic partial differential equations. This means that when the problem is studied per day and the length is measured in ????, the diffusion of substances inwards to and outwards from the tumour is almost immediate and the tumour is always in a diffusive equilibrium state.

2.1.2. Derivation of the Boundary Conditions

Let us consider an elementary test cylinder lying symmetrically across the tumour interface ???? with ??=??,??,??,??. Applying the mass conservation law in this cylinder, with elementary volume ???? and bases ??+??, ??-??, we follow similar arguments as in Section 2.1.1 to obtain ???+????·???+???(??)????-??-????·???-???(??)????=????G??????????,(2.10) where ??+?? and ??-?? are the nutrient concentration at the exterior and the interior region to ????, respectively. Taking the limit at which the test volume is eliminated, both its bases coincide to the corresponding elementary section on ???? and the volume integral of the right-hand side of (2.10) is eliminated, as G??/????is a continuous and bounded function. Applying then the integral mean value theorem, we obtain ??·???+??(??)-??·???-??(??)=0for???????with??=??,??,??,??.(2.11)

Finally, we assume continuity conditions of the nutrient concentration on each interface and an external nutrient supply traced on ???? as ??8??????=??8????????(??,??)for?????????.(2.12)

Similarly, we derive continuity boundary conditions for the inhibitor flow through the tumour’s interfaces and we assume that the inhibitor concentration is also continuous there. Finally, the asymptotic condition lim???8????(??)=0(2.13) is supposed to hold uniformly over all directions.

2.2. The Pressure Field Problem
2.2.1. Derivation of the Partial Differential Equations

Combining previous assumptions [6, 19], which attribute the cells movement either to the nutrient gradient or to the pressure gradient, as it is dictated by the Darcy law, we further assume that cells move towards the direction of the nutrient gradient and opposite to the direction of the pressure and inhibitor gradients. Denoting by ????(??) the velocity of the cell at the point ?? of the region O?? with ??=??,??,??,??, we have ????(??)=-?????????(??)+?????????(??)-?????????(??),(2.14) where ????, ????, and ???? are constants of proportionality characterizing the motility of the cell. Applying the divergence operator on both sides of (2.14), we obtain ?·????(??)=-?????????(??)+?????????(??)-?????????(??).(2.15)

Here we recall that both the tumour tissue and the host tissue are assumed to be incompressible fluids with constant density. The mass conservation law implies that ((??????(??,??))/????)+?·(????(??,??)????(??,??))=????(??,??), where ????(??,??) is the density of the tumour and ????(??,??) is the mass per unit volume, per unit time that is produced or lost in the region O?? for ??=??,??,??. Since the density is considered to be constant, we derive that ?·??????(??,??)=??(??,??)????(??,??)=????(2.16) in the tumour regions O?? for ??=??,??,??. We assume that ???? is not time dependent and it is invariable in every point inside O??. With similar arguments, we conclude that ?·????(??,??)=0 in O??, where no cell mass is neither produced nor lost.

Substituting the partial differential (2.6) and (2.7)-(2.8) in (2.15) for each region, we conclude to the following: ???????(??)=-??????+????????????????=????for???O??,???????(??)=-??????+????????G??????-????????????????=????for???O??with??=??,??,(2.17) while ?????(??)=0for???O??.(2.18)

2.2.2. Derivation of the Boundary Conditions

The tumour colony and its host surrounding are assumed to be incompressible fluids of different phase. So the Young-Laplace equation is considered on their interface ????, that is, ????(??)-??????(??)=??????for???????,(2.19) where ???? is a positive proportionality constant. Moreover, both the normal component and the tangential component of the pressure gradient are considered to be continuous across ???? in order to secure the tumour’s boundary compactness. Between the tumour’s inner interfaces we consider continuity conditions for both the pressure field and its normal derivatives, as the tumour’s distinct regions are considered to correspond to cells at different stage of their biological cycle but not to different fluids, characterized by different physical parameters. The tangential derivatives of the pressure field on the exterior surface ???? are also considered to be continuous, so that it maintains its compactness. Finally, at ????, the pressure field’s trace is ??8??????=??8??????h(??,??)for?????????.(2.20)

2.3. The Evolution Equation and Relevant Topics

The aim of this work is to determine the evolution of the tumour’s exterior boundary ???? under the assumption that this boundary evolves normally to itself. Considering that ????(????)=(??????)/???? with ??=??,??,??,??, then directly from (2.14), we easily get the following relationship: ??·??????????=-??????·???????????+??????·???????????-??????·???????????(2.21) or ??????=???????????-??????????????+??????????????-?????????????.??(2.22)

We note here that relationship (2.22) is an ordinary differential equation with respect to the function ????(??). The uniqueness of its solution is secured by the initial condition ????(0)=????, where ???? is the initial radius of the ???? boundary.

The right-hand side of (2.22) depends, as it will be shown in what follows, on the time-dependent boundaries ????(??), ????(??), ????(??), and ????(??). Hence, (2.22) is solvable under constraints, which interrelate these boundaries and secure that (2.22) is dependent only on ????(??). These constraints are provided by the critical values of the nutrient and inhibitor concentrations. In particular, the critical nutrient value ??*?? determines if a cell dies out of starvation or not, so this value is met on the surface ????, that is, ??????????=??????????=??*??.(2.23)

Also, the nutrient value ??*?? and the inhibitor value ??*?? determine if a cell proliferates or not, so these critical values are met on surface ????, that is, ??????????=??????????=??*??,??(2.24)????????=??????????=??*??.(2.25)

3. Model Analysis

In what follows the boundary value problems for the nutrient concentration field, for the inhibitor concentration, and for the pressure field are solved, using standard spectral analysis of the Laplace operator. The inhomogeneity imposed by the nutrient supply is reflected upon the conditions (2.23) and (2.24), which are considered to hold either globally on the whole surfaces ???? and ????, respectively, or as an average. Both aspects lead to the same evolution equation but allow different distribution of the nutrient supply and of the pressure on ????. The latter approach is illustrated by an example in Section 4.1, where a special inhomogeneity of the nutrient supply is assumed.

3.1. The Nutrient Concentration Problem

With respect to Section 2.1, the boundary value problem that determines the nutrient concentration field assumes the form ?????(??)=0for???O??with??=??,??,(3.1) while ?????G(??)=??????=????for???O??with??=??,??,(3.2) where the boundary conditions to be satisfied are ????(??)=????(??)for???????,????????????(??)-????????(??)=0for???????,(3.3) while ????(??)=????(??)for???????,????????????(??)-????????(??)=0for???????,????(??)=????(??)for???????,????????????(??)-????????(??)=0for???????,(3.4) whilst at the outer boundary we obtain ??8??????=??8??????????????=??8????????(??,??)for?????????.(3.5)

Before we proceed with the solution of the problem, it is worth mentioning the eigenfunctions we will use, which are the surface spherical harmonics, defined as ??????,??(??,??)=???????(cos??)cos????,??=??sin????,??=??,(3.6) which satisfy the following orthogonality property: ???2??????,??(??,??)????',??'??'1(??,??)????(??,??)=????4??(2??+1)(??+??)!??(??-??)!????'??????'??????',where????=?1,??=02,??=1.(3.7)

Moreover, for notation convenience, we will replace the triple summation symbol by the following ?+8??=0?????=0???=??,????=??,??,??? in all the expansions.

Within this aim, we expand the nutrient concentration in each domain, in spherical harmonics, demanding the field to be continuous at the center of the spherical core. In that way we, obtain ????(???)=??,??,????(??)??,????,????????????,??(??,??)for???O??,??????(??)=??6??2+???,??,?????(??)??,????,??????+??(??)??,????,????-(??+1)???????,??(??,??)for???O??,??????(??)=??6??2+???,??,?????(??)??,????,??????+??(??)??,????,????-(??+1)???????,??(??,??)for???O??,?????(??)=??,??,?????(??)??,????,??????+??(??)??,????,????-(??+1)???????,??(??,??)for???O??.(3.8)

In addition, the spherical expansion of the nutrient supply field ??8(????) on ???? reads as ??8??????=??8????????(??,??)=??8?????????,??,????????,????????,??(??,??),(3.9) where ??????,??=((2??+1)/4??)((??-??)!/(??+??)!)?????02?????0??(??,??)??????,??(??,??)sin??????????. Applying boundary conditions (3.3)–(3.5) into expansions (3.8) with (3.9), one has to deal with a seventh-order linear algebraic system with respect to the seven sequences of unknown coefficients. Long but straightforward calculations lead to the nutrient concentration field, where ????(??,??,??) is defined for ???O?? with ??=??,??,??,??, that is, ??????(??,??,??)=??-????2??2???21-3?????????+????2??2???21-3?????????-????2??2???21-3?????????+??8?????????,??,????????,??????????????????,????(??,??),(3.10)????(??,??,??)=??-????2??2???21-3?????????-????2??2???21-3?????????+????3??3???1??-1?????+????6??2+??8?????????,??,????????,??????????????????,????(??,??),(3.11)????(??,??,??)=??6??2-????2??2???21-3?????????+?????3??3??+????-????3??3??1????-1?????+??8?????????,??,????????,??????????????????,????(??,??),(3.12)?????(??,??,??)=??3??3??+????-????3??3??-????3??3??1????-1?????+??8?????????,??,????????,??????????????????,??(??,??).(3.13)

Substituting expression (3.10) into condition (2.23), we perform trivial calculations to obtain the following equation: ??*??=????-????2??2???21-3?????????-????2??2???21-3?????????+????3??3???1????-1?????+????6??2??+??8?????????,??,????????,????????????????????,??(??,??).(3.14)

Equation (3.14) can be considered in two ways, either as one that holds true pointwise on all (????,??,??)????? or as an equation that holds true at the average over the surface ????.

In the first case, we imply the orthogonality property of ??????,??(??,??) in (3.14) and we conclude that all coefficients ??????,?? for ??=1 should vanish and the nutrient supply should be defined uniformly over all directions on ????, meaning that ??8(????)=??8(????), where we have assumed for simplicity and without loss of generality that ??00,??=1.

In the second case, we take the mean value of both sides of (3.14) by integrating over ????. Due to orthogonality of ??????,??(??,??), the angular dependence will vanish and the result will read as ??*??=????-????2??2???21-3?????????-????2??2???21-3?????????+????3??3???1????-1?????+????6??2??+??8??????,(3.15) which is also the result if the first argument is considered. The difference lies in that, considering (3.14) as an average equation, the external nutrient supply is allowed to be defined nonuniformly. Though the drawback in this case lies on the perturbation implied in the locus of the points that enjoy the critical nutrient concentration ??*??, which is no longer a spherical surface.

In what follows, the treatment of either case is indifferent. Nevertheless, we will come back to the different approaches in Section 3.4, when the evolution equation is studied. Therefore, relationship (2.24) along with (3.12) gives in a similar way ??*??=????-????2??2???21-3?????????-????2??2???21-3?????????+????3??3???1????-1?????+????6??2??+??8??????.(3.16)

Both (3.15) and (3.16) should hold throughout all time of the tumour’s evolution, which provides restrictions on the evolution of all boundaries ????(??), ????(??), ????(??) and ????(??). In addition, from relations (3.15) and (3.16), we are provided with a relation depending only on ????(??) and ????(??), that is, ??*??-??*??=????3??3???1????-1?????+????6???2??-??2???.(3.17)

3.2. The Inhibitor Concentration Problem

The spherical expansions of the inhibitor concentration that solve the boundary value problem (2.7)-(2.8), with continuous boundary conditions and also with (2.13), assume the following forms: ??????(??)=??6??2+???,??,????(??)??,????,????????????,??(??,??)for???O??,??????(??)=??6??2+???,??,?????(??)??,????,??????+??(??)??,????,????-(??+1)???????,??(??,??)for???O??,??????(??)=??6??2+???,??,?????(??)??,????,??????+??(??)??,????,????-(??+1)???????,??(??,??)for???O??,?????(??)=??,??,????(??)??,????,????-(??+1)??????,??(??,??)for???R3-?O???O???O???,(3.18) where for simplicity we have used the notation -(????/????)=???? and (????/????)=???? with ??=??,??.

Following similar track of calculations as in Section 3.1, we arrive at the following expression for the inhibitor concentration field: ??????(??)=??-????2??2??+????-????2??2??-????2??2??+????6??2,??(3.19)????(??)=??-????2??2??-????2??2??+????-????3??3????+????6??2,??(3.20)????(??)=-??2??2??+?????-????3??3??+????-????3??3???1??+????6??2??,(3.21)?????(??)=??-????3??3??+????-????3??3??-????3??3???1??.(3.22)

The inhibitor concentration is, by definition, a decreasing function of ?? in R3-O??. The critical inhibitor concentration ??*?? appears on ????, and its further decrease permits cell proliferation in O??. Hence, (2.25) together with (3.20) leads to ??*??=?????2-????3???2??-????2??2??+????-????3??3??????,(3.23) which provides another restriction between the boundaries ????(??), ????(??), and ????(??), which should hold for every ??>0.

3.3. The Pressure Field

The pressure field satisfies the boundary value problem (2.17)–(2.20) completed with continuity conditions for both the field and its gradient in all inner interfaces. The corresponding spherical expansions read as ??????(??)=??6??2+???,??,????(??)??,????,????????????,??(??,??)for???O??,??????(??)=??6??2+???,??,?????(??)??,????,??????+??(??)??,????,????-(??+1)???????,??(??,??)for???O??,??????(??)=??6??2+???,??,?????(??)??,????,??????+??(??)??,????,????-(??+1)???????,??(??,??)for???O??,?????(??)=??,??,?????(??)??,????,??????+??(??)??,????,????-(??+1)???????,??(??,??)for???O??,(3.24) while the exterior pressure distribution ??8(????) on ???? is expanded in spherical harmonics as ??8??????=??8??????h(??,??)=??8?????????,??,??h????,????????,??(??,??),(3.25) where h????,??=((2??+1)/4??)((??-??)!/(??+??)!))?????02?????0h(??,??)??????,??(??,??)sin??????????. The solution of the above boundary value problem is obtained in the form ??????(??)=??-????2??2???21-3?????????+????-????2??2???21-3?????????-????2??2???21-3?????????+????????+????6??2+??8?????????,??,??h????,??????????????????,??(????,??),(3.26)????(??)=??-????2??2???21-3?????????-????2??2???21-3?????????+????????+????-????3??3???1??-1?????+????6??2+??8?????????,??,??h????,??????????????????,????(??,??),(3.27)????(??)=-??2??2???21-3?????????+????????+?????-????3??3??+????-????3??3??1????-1?????+????6??2+??8?????????,??,??h????,??????????????????,????(??,??),(3.28)?????(??)=??-????3??3??+????-????3??3??-????3??3??1????-1?????+??8?????????,??,??h????,??????????????????,??(??,??).(3.29)

3.4. The Evolution Equation

In the sequel, we collect all the results for the concentration fields and for the pressure field, as traced on ????, and we substitute them into the differential equation (2.22). Here, we have to recall the two different approaches that are possible for the nutrient concentration field, mentioned in Section 3.1.

3.4.1. The Homogeneous Approach

Following the demand that the nutrient trace on the exterior boundary ???? is homogeneous, the nutrient concentration (3.12) in the proliferating region is ????(??,??,??)=??8??????+????6??2-????2??2???21-3?????????+?????3??3??+????-????3??3??1????-1?????.(3.30)

Substituting (3.30), (3.21), and (3.28) in (2.22) and performing algebraic manipulations the following evolution equation is obtained: ??????????=???????????????3?????????-????3+????????-????3-????????3?+??????????3?????????-????3+????????-????3-????????-????3?-????????3-????????3+????????3?-??????8?????????????,??,??h????,??????????????????????,??(??,??).(3.31)

We note here that the whole model’s treatment is based on the assumption that as the tumour evolves, all its boundaries remain members of the concentric spherical coordinate system, so that all the Laplace and the Poisson partial differential equations can be treated analytically through the separation of variables method. This implies that ????=????(??), and consequently its derivative cannot depend on the direction ??=(??,??). As a consequence, only the zeroth-order term of the spherical expansion can be included in (3.31), a fact that eliminates the effect of any pressure that is imposed to the host tissue by the exterior medium on the tumour’s evolution. Moreover, the host tissue’s boundary ????=????(??) itself does not seem to be involved by the tumour’s evolution. Therefore, the evolution equation is written as follows: ??????????=???????????????3?????????-????3+????????-????3-????????3?+??????????3?????????-????3+????????-????3-????????-????3?-????????3-????????3+????????3?,(3.32) where the physical implementation was secured mathematically by imposing h????,??=0 for ??=1, ??=0,1,,??, and ??=??,?? into (3.31) concluding to (3.32). This equation describes the transitory phase until it reaches its steady state. From (3.32), one can obtain the radii of the tumour’s steady-state boundaries by setting the right-hand side equal to zero.

3.4.2. The Nonhomogeneous Approach

By considering (3.14) as an average condition, the nutrient concentration in the proliferation shell reads as ??????(??,??,??)=??6??2-????2??2???21-3?????????+?????3??3??+????-????3??3??1????-1?????+??8?????????,??,????????,??????????????????,??(??,??).(3.33)

Then, the corresponding evolution equation for ??????/???? becomes ??????????=???????????????3?????????-????3+????????-????3-????????3?+??????????3?????????-????3+????????-????3-????????-????3?-????????3-????????3+????????3?-??????8?????????????,??,??h????,??????????????????????,??(??,??)+??????8?????????????,??,????????,??????????????????????,??(??,??),(3.34) or ??????????=???????????????3?????????-????3+????????-????3-????????3?+??????????3?????????-????3+????????-????3-????????-????3?-????????3-????????3+????????3?+1???????,??,???????????????????????8????????????,??-??????8??????h????,?????????,??(??,??).(3.35)

The arguments stated in Section 3.4.1 hold true in this case as well, provided that between the spherical coefficients of the nutrient exterior field and the pressure field, as they are dictated on ????, the following condition is considered: ??????,??=??????8????????????8??????h????,??for??=1,2,,??=1,2,,??and??=??,??,(3.36) and hence ??8??????=??8??????-??????????8??????h00,??+??????????8??????(3.37) or ???????8??????-??8???????=???????8??????-??8??????h00,???.(3.38)

Consequently, the tumour’s boundaries evolve as concentric spheres and the evolution (3.35), provided (3.36), coincides with (3.32). From there after, it is obvious that the nonhomogeneous and the homogeneous approaches are treated in the same way.

We stress that the restrictions on the type of inhomogeneity exhibited by the exterior pressure are implied by the assumption (2.22) for the evolution equation, in the model presented in this work. Nevertheless, following a different approach on the evolution condition, one could investigate further the possibility of imposing an inhomogeneous pressure trace upon ????, which is under our current investigation. For this reason, we do not imply this restriction on the results of Section 3.3, since they hold true up to the point that they are imposed properly in Section 3.4 and so after.

Equation (3.32) is a nonlinear ordinary differential equation, which includes the rates (????/????) and (????/????), accompanied by ????, which are all unknown functions of time. In order to make (3.32) solvable, we need the expressions that connect all the unknown rates with function ????=????(??). This is accomplished via expressions (3.15)–(3.17) and (3.23). In particular, (3.17) can be written as ??*??-??*????2??=????3??????????3?????????-?????????+????6???????????2-??????????2?,(3.39) or using the following simple symbols for the rates under consideration, those are (????/????)=??=??(??), (????/????)=??=??(??), and (????/????)=??=??(??), where by definition 0<??,??,??<1, expression (3.39) yields -????3??3+????2??2????-??6??3-??*??-??*????2????=0.(3.40)

The third-degree nonlinear algebraic equation (3.40) provides ?? as a function of ?? and ????. Similarly, (3.23) assumes the form ????-????3??3+?????2-????3???3-?????2+??*????2?????=0.(3.41)

Equations (3.40) and (3.41) form a nonlinear algebraic system, which when it is solved, provide the rates (????/????)=??(??) and (????/????)=??(??) with respect to ????(??) and the differential (3.32) is solvable, under the initial condition ????(0)=????, which corresponds to the first observed tumour radius.

Continuing, the rate (????/????) is also connected to ????=????(??) via (3.15), or equivalently (3.16); it is worthwhile to specify this connection, as it provides the evolution of the host boundary and also it appears in the nutrient field, the inhibitor field, and the pressure field provided in Sections 3.1, 3.2, and 3.3, respectively. In particular, using (3.15), we arrive at ??*??-??8????????2??=????-????2??????????2?21-3?????????????????-????2?21-3?????????+????3??????????3?????????-?????????+????6??????????2(3.42) or, in terms of the ??,??,?? symbols, ?-????3??3-????-????3??3+????3?????=-??2??2-????-????2??2+????2+??*??-??8????????2??,(3.43) and if one substitutes the (??,??) solution of the system (3.40) and (3.41), then (3.43) provides the rate (????/????)=??(??) as a function of ????(??).

4. Special Cases

In this section, we implement the results obtained in Section 3 in three special cases that correspond to different physical approaches. At the end, we reduce the results to recover the corresponding existing results that appear in the relative literature, mainly in [5, 6].

4.1. Inhomogeneous Data Imposed on ????

For the purpose of illuminating the inhomogeneous case, we suppose that the pressure’s trace on ???? implies an inhomogeneity in the form ??8??????=??8???????1±??????????????sin2???for?????????,(4.1) where ???? is a positive proportionality constant. Expression (4.1) in terms of spherical harmonics assumes the form ??8??????=??8???????1±2????3??????????±2????3????????????2?(cos??),(4.2) which in terms of (3.25) corresponds to h00,??=?1±2????3?????????,h??20,??=52?4??02?????0?±2????3????????????2???(cos??)20,??(??,??)sin??????????=±2????3??????????,(4.3) and h20,??=0, h2??,??=0 for ??=1,2 and ??=??,??, while h????,??=0 for every ??=1,3,, ??=0,1,,?? and ??=??,??. Then, the pressure field (3.26)–(3.29) enjoys a corresponding modification. In order for the evolution of the tumour’s exterior boundary to maintain angular independence, a solution of the model exists, when the nutrient supply is also nonhomogeneous and its spherical coefficients are defined by (3.36). Then, ??8??????=??8??????+??8?????????,??,????=1??????8????????????8??????h????,????????,??(??,??)(4.4) or, for the case under consideration, ??8??????=??8??????+??????????8???????±2????3????????????2?(cos??).(4.5)

The evolution equation is not affected by this inhomogeneity, and it is given by (3.32).

4.2. Homogeneous Pressure Distribution Implying the Young-Laplace Law on ????

In this Section we consider the host nutritive environment to lie within an infinite medium, considered to be an incompressible fluid characterized by physical parameters different from those of both the host tissue and the tumour. Let us suppose that this medium imposes in the host tissue a uniform pressure ??8. Then, in the interface ????, the Young-Laplace equation should hold, which yields the boundary condition ??????????-??8=????????for?????????,(4.6) where ???? is a positive proportionality constant. Condition (4.6) corresponds to the results of the problem under investigation, if the trace of the pressure field on ???? is interpreted as ??????????=??8??????=??8+????????.(4.7)

In terms of expression (2.20), which yields ??8(????)=??8+(????/????), we have h????,??=0 for ??=1, ??=0,1,,??, and ??=??,??, while h00,??=1. Then, the pressure field (3.26)–(3.29) would be affected as follows: ????(??)=??8+????????+????????+????-????2??2???21-3?????????+????-????2??2???21-3?????????-????2??2???21-3?????????+????6??2,????(??)=??8+????????+????????+????-????2??2???21-3?????????-????2??2???21-3?????????+????-????3??3???1??-1?????+????6??2,????(??)=??8+????????+????????-????2??2???21-3?????????+?????-????3??3??+????-????3??3??1????-1?????+????6??2,????(??)=??8+????????+?????-????3??3??+????-????3??3??-????3??3??1????-1?????.(4.8)

Though, the exterior pressure should not affect the evolution equation at all, as it is explained in Section 3.4.

4.3. Homogeneous Nutrient Concentration Implying the Gibbs-Thomson Relation on ????

Following a related idea presented in [11], one can assume the Gibbs-Thomson relation for connecting the nutrient concentration trace at the point ????????? with the corresponding curvature of ????, that is, ??????????-??8=????,??????for?????????,(4.9) where ????,?? is a positive proportionality constant. This jump condition is justified if we consider two facts. Firstly, as in Section 4.1, we consider that the host nutritive environment lies within an infinite medium of constant nutrient concentration ??8. Secondly, we assume that the noncancerous cells at the boundary ???? of the host organ consume nutrients in order to provide the energy necessary for maintaining the bonds between them and the solidity of the whole structure and this energy is curvature related as dictated by the Gibbs-Thomson relation (4.9). This case corresponds to the results obtained in the present paper if we imply ??8??????=??8??????=??8+????,??????(4.10) inside (3.9).

The nutrient field in this case is immediately derived from relationships (3.10)–(3.13), by substitution of (4.10), in the form ????(??)=??8+????,??????+????-????2??2???21-3?????????+????2??2???21-3?????????-????2??2???21-3?????????,??????(??)=??-????2??2???21-3?????????-????2??2???21-3?????????+????3??3???1??-1?????+????6??2+??8+????,??????,??????(??)=-??2??2???21-3?????????+?????3??3??+????-????3??3??1????-1?????+????6??2+??8+????,??????,???????(??)=??3??3??+????-????3??3??-????3??3??1????-1?????+??8+????,??????.(4.11)

The corresponding modifications upon (3.15) and (3.16), related to the nutrient critical values, are ??*??=????-????2??2???21-3?????????-????2??2???21-3?????????+????3??3???1????-1?????+????6??2??+??8+????,??????,??*??=????-????2??2???21-3?????????-????2??2???21-3?????????+????3??3???1????-1?????+????6??2??+??8+????,??????.(4.12)

Since (3.15) and (3.16) provide the restrictions that interconnect the tumour’s boundaries, the evolution equation is affected by the assumption (4.10) through its effect upon (4.12).

The above results fully recover the corresponding nutrient results obtained in [11] for the case where ????=???? and ????=????.

4.4. Recovery of Existing Results

In order to let the results obtained in Section 3 to recover existing results in the literature [5, 6], we take the limit of our results where ?????????, at the case where ????=???? and ????(????)=??8. Then, expressions (3.10)–(3.13) become ??????(??,??,??)=??2??2???21-3?????????-????6??2??+??8,????(??,??,??)=??????(??,??,??)=-??6??2??+????3??3???1??-1?????+????6??2+??8.(4.13)

Moreover, if the inhibitor is considered to vanish far away from the tumour and not on the tumour’s boundary and also if ????=????=0, then the results (3.19)–(3.22) for the inhibitor concentration coincide with the case examined in [5] and yield ??????(??)=-??2??2??+????6??2,????(??)=??????(??)=-??3??3????.(4.14)

Then, (3.15) and (3.23) reveal ??8-??*????=-??3?12???2??-??2???-??3???1????-1????,????*????=-??3??3??????,(4.15) respectively, while they provide the inner boundaries ????, ???? with respect to the exterior boundary ???? of the tumour.

The corresponding results for the pressure field when no exterior pressure is imposed on the tumor are found in [6], where the interior of the tumor is characterized by ????=????=???? and ??8(????)=0. Then, the expressions (3.26)–(3.28) meet the corresponding result in [6] and become ????(??)=????(??)=??????(??)=??6???2-??2???+????????.(4.16)

The evolution equation studied in [6] coincides with the one included in the present work with the parameter values ????=1, ????=(??/??), and ????=0. The results cannot be compared further with [5] or [6] as they have been obtained either due to the use of different evolution equation or because the nutrient satisfies different conditions on the tumour’s boundary.

5. Discussion

In the present work, we have studied a continuous model of avascular tumour growth, under two basic assumptions. Firstly, it evolves maintaining a spherical multilayer structure, lying inside a finite concentric spherical host medium. Secondly, its evolution is regulated by the diffusion of an inhomogeneous nutrient field and of a growth inhibitory factor and by the pressure that results from the cell proliferation and disintegration, as well as from an externally imposed inhomogeneous pressure from the finite surrounding medium. All the structure interfaces evolve in time, but the time scale of the space expansion of the tumour is so much different than the diffusion time scales of the substances in the tumour that the equilibrium assumption is well justified.

Hence, the model is formulated in three boundary value problems that hold true as the tumour evolves and provide the nutrient field, the internally produced inhibitor field and the pressure field throughout the tumour and the host surrounding. The model includes an assumption for the boundaries’ evolution, which is formulated as a nonlinear ordinary differential equation with respect to the tumour’s exterior boundary, and also it includes connection formulae between all the other boundaries with respect to the tumour’s exterior one.

The same formulation and analysis can probably be adapted in different models with alternative interpretations, for example, the nutrient concentration can be considered as a drug substance supplied externally that inhibits the growth, in which case the evolution equation would be respectively modified and so on.

In the present work we focus on the type of inhomogeneity that can be compatible with the specific model that we present. The inhomogeneities that we are focused on concern the data that are independent of the tumour’s development, that is, the nutrient supply and the pressure field imposed by the exterior to the host organ surrounding. It turns out that a concentric spherical multilayer development could be secured under homogeneous nutrient supply and pressure field, or by a nutrient supply and an exterior pressure that exhibit a special connection between them, expressed in (3.36). In particular, the assumption that the spherical interfaces that separate the tumour regions are characterized by the critical nutrient values, demands that the nutrient supply should be homogeneous. Moreover, the assumption that the interfaces evolve maintaining their concentric spherical shape, together with the specific assumption for the tumour’s evolution, impose, the direction independence for the pressure trace on the exterior host tissue surface. On the other hand, relaxing the first demand and allowing the interfaces to be perturbed, the nutrient supply can be inhomogeneous and consequently the pressure trace can exhibit a type of inhomogeneity that secures the angular independence of the evolution equation.

These considerations are illuminated, by implementing the general results in special cases, in Section 4. In particular, it is shown that an angular-dependent pressure field of the form (4.1) permits a concentric spherical tumour development only under a similarly nonhomogeneous nutrient distribution in the surrounding, given in (4.5). A homogeneous pressure of the form (4.7), that takes into account the exterior boundary’s curvature and is dictated by the Young-Laplace law of different fluids’ interface, is shown to have no effect on the evolution equation of the tumour’s boundary and also to be irrelevant of the form of the nutrient supply. On the other hand, a homogeneous nutrient field that takes into account the exterior spherical curvature, dictated by the Gibbs-Thomson relation, has an indirect effect on the tumour’s evolution, as it effects the connection formulae (4.12), between the tumour’s interfaces. Though, as the two latter cases refer to homogeneous imposed fields, the interrelation (3.36), in view of (3.9) and (3.25), does not exist, as it is expected. Finally, the reduction to the simple cases of both homogeneous fields in an infinite host surrounding exhibits the lack of interdependence between the imposed pressure and the nutrient supply and fully recovers the existing radial results from the literature, as mentioned in Section 4.4. The consideration of the model presented in Section 3 under shape perturbations is not examined in the present research, as the focus of this work is on studying the consistency of the spherical development with nonhomogeneous imposed data. However, the sensitivity analysis of the proposed model is within our future consideration.

Clearly, it is a drawback of the model presented, that it incorporates a rather ideal concentric spherical structure and the results refer to such an idealistic configuration or simplified assumptions. Alternative evolution approaches for the same spherical structure in avascular tumour development, and also alternative geometrical structure of the development, which is much more applicable to cancer growing in humans, are under our current investigation.