About this Journal Submit a Manuscript Table of Contents
International Journal of Differential Equations
Volume 2012 (2012), Article ID 175434, 25 pages
http://dx.doi.org/10.1155/2012/175434
Research Article

The Avascular Tumour Growth in the Presence of Inhomogeneous Physical Parameters Imposed from a Finite Spherical Nutritive Environment

1School of Science and Technology, Hellenic Open University, 262 22 Patras, Greece
2Department of Engineering Sciences, University of Patras, 265 04 Patras, Greece

Received 20 April 2012; Accepted 22 July 2012

Academic Editor: Shaher Momani

Copyright © 2012 Foteini Kariotou and Panayiotis Vafeas. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 Ω𝑝 is occupied by proliferating cells, the subsequent shell Ω𝑞 is occupied by quiescent cells, and the interior necrotic sphere Ω𝑛 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 𝜎(𝐫𝑒) 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 Ω𝑝. In the inner shell Ω𝑞, where 𝜎(𝐫) lies in the interval 0<𝜎𝑛<𝜎<𝜎𝑝, the tumour cells cannot proliferate and live in a quiescent state. Finally, in Ω𝑛, the concentration falls below 𝜎𝑛 and the tumour cells die out of starvation. Moreover, Ω𝑛 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 Ω𝑛, 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 Ω𝑝. 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 𝑃(𝐫𝑒), 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 𝐫3. Then, each distinct region in the tumour is defined as follows. The necrotic core is defined as Ω𝑛=(𝑟,𝜃,𝜙)30𝑟<𝑟𝑛,,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 Ω𝑞=(𝑟,𝜃,𝜙)3𝑟𝑛<𝑟<𝑟𝑞,Ω,0𝜃𝜋,0𝜙<2𝜋𝑝=(𝑟,𝜃,𝜙)3𝑟𝑞<𝑟<𝑟𝑝,Ω,0𝜃𝜋,0𝜙<2𝜋𝑒=(𝑟,𝜃,𝜙)3𝑟𝑝<𝑟<𝑟𝑒,,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.

175434.fig.001
Figure 1: The spherical domains of the tumour growth model.
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 Ω𝑖, 𝑖=𝑛,𝑞,𝑝,𝑒, the rate of change of the nutrient concentration in each region equals to the normal flow of nutrients 𝐅in(𝐫) inwards to Ω𝑖, minus the flow outwards, 𝐅out(𝐫), and minus the amount of nutrient that is consumed in the interior of Ω𝑖, that is, Ω𝑖𝜕𝜎𝑖(𝐫)𝜕𝑡𝑑𝑣=𝜕Ω𝑖̂(𝐧)𝐅in(𝐫)𝑑𝑠𝜕Ω𝑖̂𝐧𝐅out(𝐫)𝑑𝑠Ω𝑖Γ𝑖𝑑𝑣,(2.3) where 𝜕Ω𝑖 is the boundary of the region Ω𝑖, ̂𝐧 is the outward unit normal vector to 𝜕Ω𝑖, and Γ𝑖 is the nutrient consumption rate, considered to be constant in Ω𝑖. 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 Ω𝑖𝜕𝜎𝑖(𝐫)𝜕𝑡𝑑𝑣=𝑆+𝑖̂(𝐫)𝐅in(𝐫)𝑑𝑠𝑆𝑖̂(𝐫)𝐅out(𝐫)𝑑𝑠Ω𝑖Γ𝑖𝑑𝑣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 Ω𝑖𝜕𝜎𝑖(𝐫)𝜕𝑡𝑑𝑣=𝑆+𝑖̂𝑘𝐫𝜎𝜎𝑖(𝐫)𝑑𝑠𝑆𝑖̂𝑘𝐫𝜎𝜎𝑖(𝐫)𝑑𝑠Ω𝑖Γ𝑖𝑑𝑣.(2.5)

Applying the divergence theorem in Ω𝑖 in the right-hand side of (2.5), (2.5) becomes Ω𝑖(((𝜕𝜎𝑖(𝐫))/𝜕𝑡)𝑘𝜎Δ𝜎𝑖(𝐫)+Γ𝑖)𝑑𝑣=0, which results in 𝜕𝜎𝑖(𝐫)𝜕𝑡𝑘𝜎Δ𝜎𝑖(𝐫)+Γ𝑖=0for𝐫Ω𝑖with𝑖=𝑛,𝑞,𝑝,𝑒.(2.6)

In particular, considering the physical assumptions made in the Introduction, the consumption rate in each region is defined as Γ𝑒=Γ𝑛=0, in Ω𝑒 and Ω𝑛, respectively, where no consumption of nutrients takes place. In Ω𝑞 and Ω𝑝, 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 Γ𝑞=Γ(𝜎𝑛/𝜎)and Γ𝑝=Γ(𝜎𝑝/𝜎), respectively, where Γ is a proportional constant depending on the species of the host tissue cells, 𝜎𝑛, 𝜎𝑝 are the nutrient critical values mentioned in the introduction, and 𝜎 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𝐫Ω𝑛,𝜕𝛽𝑖(𝐫)𝜕𝑡𝑘𝛽Δ𝛽𝑖(𝐫)+𝐶𝑖=0for𝐫Ω𝑖with𝑖=𝑞,𝑝,(2.7) and considering that the inhibitor is neither consumed nor produced in the host tissue, we obtain 𝜕𝛽𝑒(𝐫)𝜕𝑡𝑘𝛽Δ𝛽𝑒(𝐫)=0for𝐫3Ω𝑛Ω𝑞Ω𝑝.(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 𝜎, we define the dimensionless variables 𝑡=(𝑡/𝑡𝑔),𝑟=(𝑟/𝑟𝑔)and 𝜎𝑖(𝐫)=𝜎𝑖(𝐫)/𝜎for 𝑖=𝑛,𝑞,𝑝,𝑒. Hence, (2.6) becomes 𝜀𝜕𝜎𝑖𝐫𝜕𝑡Δ𝐫𝜎𝑖𝐫+𝑟2𝑔𝜎𝑘𝜎Γ𝑖=0with𝑖=𝑛,𝑞,𝑝,𝑒,(2.9) where Δ𝐫=𝑟2𝑔Δ, while the constant 𝜀=𝑟2𝑔/(𝑡𝑔𝑘𝜎) becomes much smaller than one, when the scaling 𝑟𝑔=102 cm, 𝑡𝑔=1day=86400sec, and 𝑘𝜎=106cm2/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𝑔/(𝑡𝑔𝑘𝜎)=(104/(86400106))=1,12103, 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 𝑆+𝑖̂𝐫𝜎+𝑖(𝐫)𝑑𝑠𝑆𝑖̂𝐫𝜎𝑖(𝐫)𝑑𝑠=𝑉𝑖Γ𝑖𝑘𝜎𝑑𝑣,(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 Γ𝑖/𝑘𝜎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 𝜎𝐫𝑒=𝜎𝑟𝑒𝑔(𝜃,𝜙)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𝑟𝛽𝑒(𝐫)=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 Ω𝑖 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 Ω𝑖 for 𝑖=𝑛,𝑞,𝑝. Since the density is considered to be constant, we derive that 𝐯𝑖𝑀(𝐫,𝑡)=𝑖(𝐫,𝑡)𝜌𝑡(𝐫,𝑡)𝐺𝑖(2.16) in the tumour regions Ω𝑖 for 𝑖=𝑛,𝑞,𝑝. We assume that 𝐺𝑖 is not time dependent and it is invariable in every point inside Ω𝑖. With similar arguments, we conclude that 𝐯𝑒(𝐫,𝑡)=0 in Ω𝑒, 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𝐫Ω𝑛,Δ𝑃𝑖𝐺(𝐫)=𝑖𝜇𝑃+𝜇𝜎𝜇𝑃Γ𝑖𝑘𝜎𝜇𝛽𝜇𝑃𝐶𝑖𝑘𝛽𝐹𝑖for𝐫Ω𝑖with𝑖=𝑞,𝑝,(2.17) while Δ𝑃𝑒(𝐫)=0for𝐫Ω𝑒.(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 𝑃𝐫𝑒=𝑃𝑟𝑒(𝜃,𝜙)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𝐫Ω𝑖with𝑖=𝑛,𝑒,(3.1) while Δ𝜎𝑖Γ(𝐫)=𝑖𝑘𝜎𝛾𝑖for𝐫Ω𝑖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 𝜎𝐫𝑒=𝜎𝑟𝑒𝑔̂𝐫𝑒=𝜎𝑟𝑒𝑔(𝜃,𝜙)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 +𝑙=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𝐫Ω𝑛,𝜎𝑞𝛾(𝐫)=𝑞6𝑟2+𝑙,𝑚,𝑠𝑏(𝑖)𝑚,𝑠𝑙,𝜎𝑟𝑙+𝑏(𝑒)𝑚,𝑠𝑙,𝜎𝑟(𝑙+1)𝑌𝑙𝑚,𝑠(𝜃,𝜙)for𝐫Ω𝑞,𝜎𝑝𝛾(𝐫)=𝑝6𝑟2+𝑙,𝑚,𝑠𝑐(𝑖)𝑚,𝑠𝑙,𝜎𝑟𝑙+𝑐(𝑒)𝑚,𝑠𝑙,𝜎𝑟(𝑙+1)𝑌𝑙𝑚,𝑠(𝜃,𝜙)for𝐫Ω𝑝,𝜎𝑒(𝐫)=𝑙,𝑚,𝑠𝑑(𝑖)𝑚,𝑠𝑙,𝜎𝑟𝑙+𝑑(𝑒)𝑚,𝑠𝑙,𝜎𝑟(𝑙+1)𝑌𝑙𝑚,𝑠(𝜃,𝜙)for𝐫Ω𝑒.(3.8)

In addition, the spherical expansion of the nutrient supply field 𝜎(𝐫𝑒) on 𝑆𝑒 reads as 𝜎𝐫𝑒=𝜎𝑟𝑒𝑔(𝜃,𝜙)=𝜎𝑟𝑒𝑙,𝑚,𝑠𝑔𝑙𝑚,𝑠𝑌𝑙𝑚,𝑠(𝜃,𝜙),(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 𝐫Ω𝑖 with 𝑖=𝑛,𝑞,𝑝,𝑒, that is, 𝜎𝑛𝛾(𝑟,𝜃,𝜙)=𝑝𝛾𝑞2𝑟2𝑞213𝑟𝑞𝑟𝑒+𝛾𝑞2𝑟2𝑛213𝑟𝑛𝑟𝑒𝛾𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝜎𝑟𝑒𝑙,𝑚,𝑠𝑔𝑙𝑚,𝑠𝑟𝑟𝑒𝑙𝑌𝑙𝑚,𝑠𝜎(𝜃,𝜙),(3.10)𝑞𝛾(𝑟,𝜃,𝜙)=𝑝𝛾𝑞2𝑟2𝑞213𝑟𝑞𝑟𝑒𝛾𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝛾𝑞3𝑟3𝑛1𝑟1𝑟𝑒+𝛾𝑞6𝑟2+𝜎𝑟𝑒𝑙,𝑚,𝑠𝑔𝑙𝑚,𝑠𝑟𝑟𝑒𝑙𝑌𝑙𝑚,𝑠𝜎(𝜃,𝜙),(3.11)𝑝𝛾(𝑟,𝜃,𝜙)=𝑝6𝑟2𝛾𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝛾𝑞3𝑟3𝑛+𝛾𝑝𝛾𝑞3𝑟3𝑞1𝑟1𝑟𝑒+𝜎𝑟𝑒𝑙,𝑚,𝑠𝑔𝑙𝑚,𝑠𝑟𝑟𝑒𝑙𝑌𝑙𝑚,𝑠𝜎(𝜃,𝜙),(3.12)𝑒𝛾(𝑟,𝜃,𝜙)=𝑞3𝑟3𝑛+𝛾𝑝𝛾𝑞3𝑟3𝑞𝛾𝑝3𝑟3𝑝1𝑟1𝑟𝑒+𝜎𝑟𝑒𝑙,𝑚,𝑠𝑔𝑙𝑚,𝑠𝑟𝑟𝑒𝑙𝑌𝑙𝑚,𝑠(𝜃,𝜙).(3.13)

Substituting expression (3.10) into condition (2.23), we perform trivial calculations to obtain the following equation: 𝜎𝑛=𝛾𝑝𝛾𝑞2𝑟2𝑞213𝑟𝑞𝑟𝑒𝛾𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝛾𝑞3𝑟3𝑛1𝑟𝑛1𝑟𝑒+𝛾𝑞6𝑟2𝑛+𝜎𝑟𝑒𝑙,𝑚,𝑠𝑔𝑙𝑚,𝑠𝑟𝑛𝑟𝑒𝑙𝑌𝑙𝑚,𝑠(𝜃,𝜙).(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 𝜎(𝐫𝑒)=𝜎(𝑟𝑒), 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𝑞213𝑟𝑞𝑟𝑒𝛾𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝛾𝑞3𝑟3𝑛1𝑟𝑛1𝑟𝑒+𝛾𝑞6𝑟2𝑛+𝜎𝑟𝑒,(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𝑞213𝑟𝑞𝑟𝑒𝛾𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝛾𝑞3𝑟3𝑛1𝑟𝑞1𝑟𝑒+𝛾𝑞6𝑟2𝑞+𝜎𝑟𝑒.(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𝐫Ω𝑛,𝛽𝑞𝑐(𝐫)=𝑞6𝑟2+𝑙,𝑚,𝑠𝑏(𝑖)𝑚,𝑠𝑙,𝛽𝑟𝑙+𝑏(𝑒)𝑚,𝑠𝑙,𝛽𝑟(𝑙+1)𝑌𝑙𝑚,𝑠(𝜃,𝜙)for𝐫Ω𝑞,𝛽𝑝𝑐(𝐫)=𝑝6𝑟2+𝑙,𝑚,𝑠𝑐(𝑖)𝑚,𝑠𝑙,𝛽𝑟𝑙+𝑐(𝑒)𝑚,𝑠𝑙,𝛽𝑟(𝑙+1)𝑌𝑙𝑚,𝑠(𝜃,𝜙)for𝐫Ω𝑝,𝛽𝑒(𝐫)=𝑙,𝑚,𝑠𝑑(𝑒)𝑚,𝑠𝑙,𝛽𝑟(𝑙+1)𝑌𝑙𝑚,𝑠(𝜃,𝜙)for𝐫3Ω𝑝Ω𝑞Ω𝑛,(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 3Ω𝑛. The critical inhibitor concentration 𝛽𝑝 appears on 𝑆𝑞, and its further decrease permits cell proliferation in Ω𝑝. 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𝐫Ω𝑛,𝑃𝑞𝐹(𝐫)=𝑞6𝑟2+𝑙,𝑚,𝑠𝑏(𝑖)𝑚,𝑠𝑙,𝑃𝑟𝑙+𝑏(𝑒)𝑚,𝑠𝑙,𝑃𝑟(𝑙+1)𝑌𝑙𝑚,𝑠(𝜃,𝜙)for𝐫Ω𝑞,𝑃𝑝𝐹(𝐫)=𝑝6𝑟2+𝑙,𝑚,𝑠𝑐(𝑖)𝑚,𝑠𝑙,𝑃𝑟𝑙+𝑐(𝑒)𝑚,𝑠𝑙,𝑃𝑟(𝑙+1)𝑌𝑙𝑚,𝑠(𝜃,𝜙)for𝐫Ω𝑝,𝑃𝑒(𝐫)=𝑙,𝑚,𝑠𝑑(𝑖)𝑚,𝑠𝑙,𝑃𝑟𝑙+𝑑(𝑒)𝑚,𝑠𝑙,𝑃𝑟(𝑙+1)𝑌𝑙𝑚,𝑠(𝜃,𝜙)for𝐫Ω𝑒,(3.24) while the exterior pressure distribution 𝑃(𝐫𝑒) on 𝑆𝑒 is expanded in spherical harmonics as 𝑃𝐫𝑒=𝑃𝑟𝑒(𝜃,𝜙)=𝑃𝑟𝑒𝑙,𝑚,𝑠𝑙𝑚,𝑠𝑌𝑙𝑚,𝑠(𝜃,𝜙),(3.25) where 𝑙𝑚,𝑠=((2𝑙+1)/4𝜋)((𝑙𝑚)!/(𝑙+𝑚)!))𝜀𝑚02𝜋𝜋0(𝜃,𝜙)𝑌𝑙𝑚,𝑠(𝜃,𝜙)sin𝜃𝑑𝜃𝑑𝜙. The solution of the above boundary value problem is obtained in the form 𝑃𝑛𝐹(𝐫)=𝑞𝐹𝑛2𝑟2𝑛213𝑟𝑛𝑟𝑒+𝐹𝑝𝐹𝑞2𝑟2𝑞213𝑟𝑞𝑟𝑒𝐹𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝑎𝑝𝑟𝑝+𝐹𝑛6𝑟2+𝑃𝑟𝑒𝑙,𝑚,𝑠𝑙𝑚,𝑠𝑟𝑟𝑒𝑙𝑌𝑙𝑚,𝑠(𝑃𝜃,𝜙),(3.26)𝑞𝐹(𝐫)=𝑝𝐹𝑞2𝑟2𝑞213𝑟𝑞𝑟𝑒𝐹𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝑎𝑝𝑟𝑝+𝐹𝑞𝐹𝑛3𝑟3𝑛1𝑟1𝑟𝑒+𝐹𝑞6𝑟2+𝑃𝑟𝑒𝑙,𝑚,𝑠𝑙𝑚,𝑠𝑟𝑟𝑒𝑙𝑌𝑙𝑚,𝑠𝑃(𝜃,𝜙),(3.27)𝑝𝐹(𝐫)=𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝑎𝑝𝑟𝑝+𝐹𝑞𝐹𝑛3𝑟3𝑛+𝐹𝑝𝐹𝑞3𝑟3𝑞1𝑟1𝑟𝑒+𝐹𝑝6𝑟2+𝑃𝑟𝑒𝑙,𝑚,𝑠𝑙𝑚,𝑠𝑟𝑟𝑒𝑙𝑌𝑙𝑚,𝑠𝑃(𝜃,𝜙),(3.28)𝑒𝐹(𝐫)=𝑞𝐹𝑛3𝑟3𝑛+𝐹𝑝𝐹𝑞3𝑟3𝑞𝐹𝑝3𝑟3𝑝1𝑟1𝑟𝑒+𝑃𝑟𝑒𝑙,𝑚,𝑠𝑙𝑚,𝑠𝑟𝑟𝑒𝑙𝑌𝑙𝑚,𝑠(𝜃,𝜙).(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 𝜎𝑝(𝑟,𝜃,𝜙)=𝜎𝑟𝑒+𝛾𝑝6𝑟2𝛾𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝛾𝑞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𝜇𝑝𝑃𝑟𝑒𝑟𝑝𝑙,𝑚,𝑠𝑙𝑚,𝑠𝑙𝑟𝑝𝑟𝑒𝑙𝑌𝑙𝑚,𝑠(𝜃,𝜙).(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 𝑙𝑚,𝑠=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𝑝213𝑟𝑝𝑟𝑒+𝛾𝑞3𝑟3𝑛+𝛾𝑝𝛾𝑞3𝑟3𝑞1𝑟1𝑟𝑒+𝜎𝑟𝑒𝑙,𝑚,𝑠𝑔𝑙𝑚,𝑠𝑟𝑟𝑒𝑙𝑌𝑙𝑚,𝑠(𝜃,𝜙).(3.33)

Then, the corresponding evolution equation for 𝑑𝑟𝑝/𝑑𝑡 becomes 𝑑𝑟𝑝𝑑𝑡=𝑟𝑝𝑟𝑛𝑟𝑝3𝜇𝑃𝐹𝑞𝐹𝑛3+𝜇𝛽𝑐𝑞𝑝𝑛3𝜇𝜎𝛾𝑞3+𝑟𝑞𝑟𝑝3𝜇𝑃𝐹𝑝𝐹𝑞3+𝜇𝛽𝑐𝑝𝑐𝑞3𝜇𝜎𝛾𝑝𝛾𝑞3𝜇𝑃𝐹𝑝3𝜇𝛽𝑐𝑝3+𝜇𝜎𝛾𝑝3𝜇𝑝𝑃𝑟𝑒𝑟𝑝𝑙,𝑚,𝑠𝑙𝑚,𝑠𝑙𝑟𝑝𝑟𝑒𝑙𝑌𝑙𝑚,𝑠(𝜃,𝜙)+𝜇𝜎𝜎𝑟𝑒𝑟𝑝𝑙,𝑚,𝑠𝑔𝑙𝑚,𝑠𝑙𝑟𝑝𝑟𝑒𝑙𝑌𝑙𝑚,𝑠(𝜃,𝜙),(3.34) or 𝑑𝑟𝑝𝑑𝑡=𝑟𝑝𝑟𝑛𝑟𝑝3𝜇𝑃𝐹𝑞𝐹𝑛3+𝜇𝛽𝑐𝑞𝑝𝑛3𝜇𝜎𝛾𝑞3+𝑟𝑞𝑟𝑝3𝜇𝑃𝐹𝑝𝐹𝑞3+𝜇𝛽𝑐𝑝𝑐𝑞3𝜇𝜎𝛾𝑝𝛾𝑞3𝜇𝑃𝐹𝑝3𝜇𝛽𝑐𝑝3+𝜇𝜎𝛾𝑝3+1𝑟𝑝𝑙,𝑚,𝑠𝑙𝑟𝑝𝑟𝑒𝑙𝜇𝜎𝜎𝑟𝑒𝑔𝑙𝑚,𝑠𝜇𝑝𝑃𝑟𝑒𝑙𝑚,𝑠𝑌𝑙𝑚,𝑠(𝜃,𝜙).(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: 𝑔𝑙𝑚,𝑠=𝜇𝑝𝑃𝑟𝑒𝜇𝜎𝜎𝑟𝑒𝑙𝑚,𝑠for𝑙=1,2,,𝑚=1,2,,𝑙and𝑠=𝑒,𝑜,(3.36) and hence 𝜎𝐫𝑒=𝜎𝑟𝑒𝜇𝑝𝜇𝜎𝑃𝑟𝑒00,𝑒+𝜇𝑝𝜇𝜎𝑃𝐫𝑒(3.37) or 𝜇𝜎𝜎𝐫𝑒𝜎𝑟𝑒=𝜇𝑝𝑃𝐫𝑒𝑃𝑟𝑒00,𝑒.(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 𝜎𝑛𝜎𝑟𝑒𝑟2𝑝=𝛾𝑝𝛾𝑞2𝑟𝑞𝑟𝑝2213𝑟𝑞𝑟𝑝𝑟𝑝𝑟𝑒𝛾𝑝2213𝑟𝑝𝑟𝑒+𝛾𝑞3𝑟𝑛𝑟𝑝3𝑟𝑝𝑟𝑛𝑟𝑝𝑟𝑒+𝛾𝑞6𝑟𝑛𝑟𝑝2(3.42) or, in terms of the 𝑥,𝑦,𝑧 symbols, 𝛾𝑞3𝑥3𝛾𝑝𝛾𝑞3𝑦3+𝛾𝑝3𝛾𝑧=𝑞2𝑥2𝛾𝑝𝛾𝑞2𝑦2+𝛾𝑝2+𝜎𝑛𝜎𝑟𝑒𝑟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 𝑃𝐫𝑒=𝑃𝑟𝑒1±𝑎𝑒𝑟𝑝𝑟𝑒sin2𝜃for𝐫𝑒𝑆𝑒,(4.1) where 𝑎𝑒 is a positive proportionality constant. Expression (4.1) in terms of spherical harmonics assumes the form 𝑃𝐫𝑒=𝑃𝑟𝑒1±2𝑎𝑒3𝑟𝑝𝑟𝑒2𝑎𝑒3𝑟𝑝𝑟𝑒𝑃2(cos𝜃),(4.2) which in terms of (3.25) corresponds to 00,𝑒=1±2𝑎𝑒3𝑟𝑝𝑟𝑒,20,𝑒=524𝜋02𝜋𝜋02𝑎𝑒3𝑟𝑝𝑟𝑒𝑃2𝑌(cos𝜃)20,𝑒(𝜃,𝜙)sin𝜃𝑑𝜃𝑑𝜙=2𝑎𝑒3𝑟𝑝𝑟𝑒,(4.3) and 20,𝑜=0, 2𝑚,𝑠=0 for 𝑚=1,2 and 𝑠=𝑒,𝑜, while 𝑙𝑚,𝑠=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, 𝜎𝐫𝑒=𝜎𝑟𝑒+𝜎𝑟𝑒𝑙,𝑚,𝑠𝑙1𝜇𝑝𝑃𝑟𝑒𝜇𝜎𝜎𝑟𝑒𝑙𝑚,𝑠𝑌𝑙𝑚,𝑠(𝜃,𝜙)(4.4) or, for the case under consideration, 𝜎𝐫𝑒=𝜎𝑟𝑒+𝜇𝑝𝜇𝜎𝑃𝑟𝑒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 𝑃. Then, in the interface 𝑆𝑒, the Young-Laplace equation should hold, which yields the boundary condition 𝑃𝑒𝐫𝑒𝑃=𝑎𝑒𝑟𝑒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 𝑃𝑒𝐫𝑒=𝑃𝐫𝑒𝑃+𝑎𝑒𝑟𝑒.(4.7)

In terms of expression (2.20), which yields 𝑃(𝑟𝑒)=𝑃+(𝑎𝑒/𝑟𝑒), we have 𝑙𝑚,𝑠=0 for 𝑙1, 𝑚=0,1,,𝑙, and 𝑠=𝑒,𝑜, while 00,𝑒=1. Then, the pressure field (3.26)–(3.29) would be affected as follows: 𝑃𝑛(𝐫)=𝑃+𝑎𝑒𝑟𝑒+𝑎𝑝𝑟𝑝+𝐹𝑞𝐹𝑛2𝑟2𝑛213𝑟𝑛𝑟𝑒+𝐹𝑝𝐹𝑞2𝑟2𝑞213𝑟𝑞𝑟𝑒𝐹𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝐹𝑛6𝑟2,𝑃𝑞(𝐫)=𝑃+𝑎𝑒𝑟𝑒+𝑎𝑝𝑟𝑝+𝐹𝑝𝐹𝑞2𝑟2𝑞213𝑟𝑞𝑟𝑒𝐹𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝐹𝑞𝐹𝑛3𝑟3𝑛1𝑟1𝑟𝑒+𝐹𝑞6𝑟2,𝑃𝑝(𝐫)=𝑃+𝑎𝑒𝑟𝑒+𝑎𝑝𝑟𝑝𝐹𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝐹𝑞𝐹𝑛3𝑟3𝑛+𝐹𝑝𝐹𝑞3𝑟3𝑞1𝑟1𝑟𝑒+𝐹𝑝6𝑟2,𝑃𝑒(𝐫)=𝑃+𝑎𝑒𝑟𝑒+𝐹𝑞𝐹𝑛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, 𝜎𝑒𝐫𝑒𝜎=𝑎𝑒,𝜎𝑟𝑒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 𝜎. 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 𝜎𝐫𝑒=𝜎𝑟𝑒𝜎+𝑎𝑒,𝜎𝑟𝑒(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 𝜎𝑛(𝐫)=𝜎+𝑎𝑒,𝜎𝑟𝑒+𝛾𝑝𝛾𝑞2𝑟2𝑞213𝑟𝑞𝑟𝑒+𝛾𝑞2𝑟2𝑛213𝑟𝑛𝑟𝑒𝛾𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒,𝜎𝑞𝛾(𝐫)=𝑝𝛾𝑞2𝑟2𝑞213𝑟𝑞𝑟𝑒𝛾𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝛾𝑞3𝑟3𝑛1𝑟1𝑟𝑒+𝛾𝑞6𝑟2+𝜎+𝑎𝑒,𝜎𝑟𝑒,𝜎𝑝𝛾(𝐫)=𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝛾𝑞3𝑟3𝑛+𝛾𝑝𝛾𝑞3𝑟3𝑞1𝑟1𝑟𝑒+𝛾𝑝6𝑟2+𝜎+𝑎𝑒,𝜎𝑟𝑒,𝜎𝑒𝛾(𝐫)=𝑞3𝑟3𝑛+𝛾𝑝𝛾𝑞3𝑟3𝑞𝛾𝑝3𝑟3𝑝1𝑟1𝑟𝑒+𝜎+𝑎𝑒,𝜎𝑟𝑒.(4.11)

The corresponding modifications upon (3.15) and (3.16), related to the nutrient critical values, are 𝜎𝑛=𝛾𝑝𝛾𝑞2𝑟2𝑞213𝑟𝑞𝑟𝑒𝛾𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝛾𝑞3𝑟3𝑛1𝑟𝑛1𝑟𝑒+𝛾𝑞6𝑟2𝑛+𝜎+𝑎𝑒,𝜎𝑟𝑒,𝜎𝑝=𝛾𝑝𝛾𝑞2𝑟2𝑞213𝑟𝑞𝑟𝑒𝛾𝑝2𝑟2𝑝213𝑟𝑝𝑟𝑒+𝛾𝑞3𝑟3𝑛1𝑟𝑞1𝑟𝑒+𝛾𝑞6𝑟2𝑞+𝜎+𝑎𝑒,𝜎𝑟𝑒.(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 𝜎𝑒(𝐫𝑒)=𝜎. Then, expressions (3.10)–(3.13) become 𝜎𝑛𝛾(𝑟,𝜃,𝜙)=𝑞2𝑟2𝑛213𝑟𝑛𝑟𝑝𝛾𝑝6𝑟2𝑝+𝜎,𝜎𝑞(𝑟,𝜃,𝜙)=𝜎𝑝𝛾(𝑟,𝜃,𝜙)=𝑝6𝑟2𝑝+𝛾𝑝3𝑟3𝑛1𝑟1𝑟𝑝+𝛾𝑝6𝑟2+𝜎.(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 𝜎𝜎𝑛𝛾=𝑝312𝑟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 𝑃(𝐫𝑒)=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.

References

  1. T. Roose, S. J. Chapman, and P. K. Maini, “Mathematical models of avascular tumor growth,” SIAM Review, vol. 49, no. 2, pp. 179–208, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  2. R. P. Araujo and D. L. S. McElwain, “A history of the study of solid tumour growth: the contribution of mathematical modelling,” Bulletin of Mathematical Biology, vol. 66, no. 5, pp. 1039–1091, 2004. View at Publisher · View at Google Scholar · View at Scopus
  3. D. S. Jones and B. D. Sleeman, “Mathematical modeling of avascular and vascular tumor growth,” in Advanced Topics in Scattering and Biomedical Engineering, pp. 305–331, World Scientific, Hackensack, NJ, USA, 2008. View at Publisher · View at Google Scholar
  4. A. C. Burton, “Rate of growth of solid tumours as a problem of diffusion,” Growth, Development and Aging, vol. 30, no. 2, pp. 157–176, 1966. View at Scopus
  5. H. P. Greenspan, “Models for the growth of a solid tumor by diffusion,” Studies in Applied Mathematics, vol. 52, pp. 317–340, 1972. View at Zentralblatt MATH
  6. H. P. Greenspan, “On the growth and stability of cell cultures and solid tumors,” Journal of Theoretical Biology, vol. 56, no. 1, pp. 229–242, 1976. View at Publisher · View at Google Scholar · View at Scopus
  7. J. A. Adam, “A simplified mathematical model of tumor growth,” Mathematical Biosciences, vol. 81, no. 2, pp. 229–244, 1986. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  8. J. A. Adam, “A mathematical model of tumor growth. II. Effects of geometry and spatial nonuniformity on stability,” Mathematical Biosciences, vol. 86, no. 2, pp. 183–211, 1987. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  9. H. M. Byrne and M. A. J. Chaplain, “Growth of necrotic tumors in the presence and absence of inhibitors,” Mathematical Biosciences, vol. 135, no. 2, pp. 187–216, 1996. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  10. H. M. Byrne and M. A. J. Chaplain, “Modelling the role of cell-cell adhesion in the growth and development of carcinomas,” Mathematical and Computer Modelling, vol. 24, no. 12, pp. 1–17, 1996. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  11. H. M. Byrne and M. A. J. Chaplain, “Free boundary value problems associated with the growth and development of multicellular spheroids,” European Journal of Applied Mathematics, vol. 8, no. 6, pp. 639–658, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  12. H. M. Byrne and S. A. Gourley, “The role of growth factors in avascular tumour growth,” Mathematical and Computer Modelling, vol. 26, no. 4, pp. 35–55, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  13. J. P. Ward and J. R. King, “Mathematical modelling of avascular-tumour growth II: modelling growth saturation,” IMA Journal of Mathemathics Applied in Medicine and Biology, vol. 16, no. 2, pp. 171–211, 1999. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  14. J. P. Ward and J. R. King, “Mathematical modelling of the effects of mitotic inhibitors on avascular tumour growth,” Journal of Theoretical Medicine, vol. 1, no. 4, pp. 287–311, 1999. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  15. D. L. S. McElwain and P. J. Ponzo, “A model for the growth of a solid tumor with non-uniform oxygen consumption,” Mathematical Biosciences, vol. 35, no. 3-4, pp. 267–279, 1977. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  16. D. L. S. McElwain and L. E. Morris, “Apoptosis as a volume loss mechanism in mathematical models of solid tumor growth,” Mathematical Biosciences, vol. 39, no. 1-2, pp. 147–157, 1978. View at Publisher · View at Google Scholar · View at Scopus
  17. G. Helmlinger, P. A. Netti, H. C. Lichtenbeld, R. J. Melder, and R. K. Jain, “Solid stress inhibits the growth of multicellular tumor spheroids,” Nature Biotechnology, vol. 15, no. 8, pp. 778–783, 1997. View at Publisher · View at Google Scholar · View at Scopus
  18. J. L. Gevertz, G. T. Gillies, and S. Torquato, “Simulating tumor growth in confined heterogeneous environments,” Physical Biology, vol. 5, no. 3, Article ID 036010, 2008. View at Publisher · View at Google Scholar · View at Scopus
  19. K. A. Landman and C. P. Please, “Tumour dynamics and necrosis: surface tension and stability,” IMA Journal of Mathemathics Applied in Medicine and Biology, vol. 18, no. 2, pp. 131–158, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  20. G. J. Pettet, C. P. Please, M. J. Tindall, and D. L. S. McElwain, “The migration of cells in multicell tumor spheroids,” Bulletin of Mathematical Biology, vol. 63, no. 2, pp. 231–257, 2001. View at Publisher · View at Google Scholar · View at Scopus