New formulation of the Gompertz equation to describe the kinetics of untreated tumors.

BACKGROUND: Different equations have been used to describe and understand the growth kinetics of undisturbed malignant solid tumors. The aim of this paper is to propose a new formulation of the Gompertz equation in terms of different parameters of a malignant tumor: the intrinsic growth rate, the deceleration factor, the apoptosis rate, the number of cells corresponding to the tumor latency time, and the fractal dimensions of the tumor and its contour.
METHODS: Furthermore, different formulations of the Gompertz equation are used to fit experimental data of the Ehrlich and fibrosarcoma Sa-37 tumors that grow in male BALB/c/Cenp mice. The parameters of each equation are obtained from these fittings.
RESULTS: The new formulation of the Gompertz equation reveals that the initial number of cancerous cells in the conventional Gompertz equation is not a constant but a variable that depends nonlinearly on time and the tumor deceleration factor. In turn, this deceleration factor depends on the apoptosis rate of tumor cells and the fractal dimensions of the tumor and its irregular contour.
CONCLUSIONS: It is concluded that this new formulation has two parameters that are directly estimated from the experiment, describes well the growth kinetics of unperturbed Ehrlich and fibrosarcoma Sa-37 tumors, and confirms the fractal origin of the Gompertz formulation and the fractal property of tumors.

Introduction

One of the most interesting problems of current oncology is the understanding of the growth kinetics of a malignant tumor, named TGK (TGK), which follows a sigmoidal law. The TGK analysis is equally made by means of graphs of the number of cancer cells (n) versus time t, named n(t); tumor volume (V) versus t, named V(t); and/or the tumor mass (m) versus t, named m(t). This is due to the close relationship between these three physical quantities. Additionally, the sigmoidal form of TGK has been described by different equations, such as Gompertz, Logistics, Bertalanffy-Richards, Kolmogorov-Johnson-Mehl-Avrami modified, being the Gompertz equation (GE) the most used [1-3].

Izquierdo-Kulich et al. [4] report the fractal origin of GE (see appendix A). This fractal origin has also been reported in [5-8] but in terms only of the fractal dimension D. Here, we have considered the one in [4] because it also takes into account the fractal structure of the boundary of the tumor.

In the different formulations of the GE [1-3] and in the experiment [9, 10] the starting point of TGK is considered when the initial number of tumor cells (n0) and the initial tumor volume (V0) satisfy the conditions n (t = 0) = n0 and V (t = 0) = V0, respectively. In preclinical studies, the researcher chooses n0/V0 depending on the purpose of the investigation. The time that elapses from the inoculation of the tumor cells in the host until the tumor reaches n0/V0 is named t0 [1, 3, 9]. Nevertheless, in clinics, n0/V0 corresponds to the tumor detected for the first time by the doctor by means of clinical and/or imaging methods. For this case, t0 is the time that elapses from the tumor formation in the organism (via chemical, biological and/or physical carcinogens) [10], until its detection for the first time. This supposes n0 ≥ nmed, where nmed is the minimum number of quantifiable cancer cells contained in the smallest measurable tumor volume, named Vmed (V0 ≥ Vmed). The post-inoculation time that elapses until the tumor reaches nmed/Vmed is named tmed (t0 ≥ tmed) [3].

In [4], it is considered the Gompertz equation given in Eq (1) (named GE1)

According the considerations in the previous paragraph, GE1 has two limitations: 1) n0 = 1, which means that the tumor has only one cell when it reaches V0, in contradiction with the experiment [9, 10]. 2) The maximum capacity of the tumor (n∞) depends only on α and β and not on n0 (n(t) = n∞ = e when t → ∞). From the mathematical point of view, n∞ is the upper asymptote of TGK. Nevertheless, in the preclinical, the condition t → ∞ is the post-inoculation time that elapses until the tumor reaches a certain volume, for which animals are sacrificed for ethical reasons [1]. In clinics, this condition means the time that elapses from the tumor formation in the organism until the patient dies.

Each undisturbed solid tumor histological variety, that grows in a type of syngeneic host to it, has its own natural history (only sigmoidal law), which does not depend on the selection of n0/V0, as observed in [3, 10–12]. In the experiment, once the researcher fixes n0/V0, t0 can be estimated a priori when the tumor latency time is known, named tobs (tobs < t0), which is the post-inoculation time that elapses until that the tumor is observed for the first time. In this case, the tumor is observable and palpable but not measurable. However, its size, named Vobs (V(t = tobs) = Vobs), is estimated following the methodology reported in [1, 3]. When the tumor reaches Vobs, it contains a number of cells, named nobs (n(t = tobs) = nobs).

The interest of including nobs/Vobs (nobs/Vobs < nmed/Vmed ≤ n0/V0) in GE is because an important part of vital cycle of a solid tumor occur before it is clinically detected (Vmed), as reported in [1, 3, 10]. Furthermore, a high cellular viability (≥ 95%) and a correct inoculation of the initial concentration of tumor cells (co) are guaranteed, tobs can be known a priori for a tumor histological variety that grows in a certain type of syngeneic host to it [3, 9–11].

As far as we reviewed, few experimental works report the analysis of TGK from Vobs [1, 3] and none of equations used to describe TGK includes nobs/Vobs. In addition, in the literature a relationship of α and β in terms of Df, df and nobs/Vobs has not been reported in the literature. Therefore, the aim of this paper is to propose a new formulation of the GE that includes nobs/Vobs, n0/V0, α, β, and to study the relation of these parameters with the fractal dimensions Df and df. The validity of this new mathematical formulation and the estimation of its parameters are determined from volumes of the Ehrlich and fibrosarcoma Sa-37 tumors that grow in BALB/c/Cenp mice, previously reported in [9]. Furthermore, the graphs of α versus df and β versus df/Df for different values of u2 (the constant of the velocity of apoptosis) and nobs are shown.

Methods

Conventional Gompertz equation

Eq (2), named GE2, is the conventional GE and the most used when TGK starts at n0/V0, given by

According to GE2, n∞ depends on n0, α and β (n(t) = n∞ = n0 e when t → ∞) and results from solving the ordinary differential Eq (3) with its initial condition, given by

GE2 suggests that n0 (constant in time) has to be included in Eq (A2). Tjørve and Tjørve [2] report that n0 acts as a parameter of shape (n∞ changes with n0) or location (n∞ remains constant).

Inclusion of n0 in Eq (A2)

In this topic was followed the methodology exposed in [4] and the initial number of tumor cells at t = 0, named n00, was included in Eq (A2), resulting the following problem

The exact solution of Eq (3) was given by with

Two inconsistencies were found in [4]: 1) the coefficient 1.5 in the parameter α of Eq (A3) was not correct but 2/3, as in Eq (6). 2) Different types of experimental tumors with the same values of df and Df had different values of α/β (we refer to the reader see Table 1 of [4]), in contrast to Eq (A3).

Table 1

Parameters of the models for the Ehrlich tumor.

Parameters	Different formulations of Gompertz equations
	GE1	GE2	GE5	GE8
α (days-1)	0.160±0.005	0.466±0.012	0.285±0.004	0.719±0.067
β (days-1)	0.122±0.007	0.261±0.007	0.261±0.007	0.261±0.007
Vobs(α,β) (cm3)	-	-	-	0.190±0.063
u2 (days-1)	0.263±0.066	0.633±0.141	0.391±0.055	0.687±0.131
df	0.720±0.061	0.768±0.056	0.764±0.032	0.611±0.052
Df	1.467±0.410	1.404±0.346	1.583±0.836	1.023±0.192
Vobs(u2,df,Df) (cm3)	-	-	-	0.190±0.041
αc (days-1)	0.163±0.003	0.471±0.009	0.286±0.005	0.724±0.055
βc (days-1)	0.134±0.104	0.287±0.005	0.275±0.009	0.261±0.007
SE	0.215±0.006	0.884±0.021	0.088±0.021	0.089±0.021
PRESS	1.313±0.154	0.015±0.012	0.015±0.012	0.016±0.012
MPRESS	1.128±0.144	0.015±0.012	0.015±0.012	0.016±0.012
r2	0.990±0.006	0.998±0.009	0.998±0.009	0.998±0.001
ra2	0.990±0.006	0.998±0.009	0.998±0.009	0.998±0.001
RMSE (cm3)	0.214±0.006	0.088±0.021	0.087±0.021	0.088±0.022
Dmax (cm3)	0.501±0.013	0.194±0.050	0.194±0.050	0.195±0.050
eα	0.042±0.015	0.073±0.030	0.053±0.021	0.095±0.047
eβ	0.040±0.018	0.046±0.019	0.048±0.022	0.047±0.020
eVobs(α,β)	-	-	-	0.033±0.009
eu2	0.046±0.007	0.052±0.023	0.051±0.013	0.082±0.025
edf	0.071±0.011	0.072±0.019	0.070±0.021	0.073±0.020
eDf	0.325±0.075	0.415±0.068	0.761±0.108	0.054±0.014
eVobs(u2,df,Df)	-	-	-	0.032±0.008

Means ± standard deviation of parameters of the Ehrlich tumor and criteria for model assessment obtained for different formulations of Gompertz equations.

Eq (5), named GE5, agrees with GE2 when . In addition, the parameters n00 and n0 coincided exactly at t = 0. The constant parameter n00 (n00 ≥ nmed) constituted the starting point of TGK for GE5 and reached for t = t0. Therefore, it was convenient to differentiate n0 and n00 to compare GE2 and GE5 in order to avoid confusion in the interpretation of these two parameters. GE5 revealed that n∞ depends only on α and β and not on n00 (n(t) = n∞ = e for t → ∞).

Inclusion of nobs in GE

Eq (3) was rewritten as where n000 was the number of tumor cells that the researcher selected at t = t0. The analytical solution of Eq (7) was given by

Eq (8), named GE8, agreed with GE5 at t = 0 (for all nobs) and when nobs = 1 (for all t). The GE8 coincided with the GE2 at t = 0 (for all nobs) and when . The parameter nobs (nobs < nmed ≤ n000) was the starting point of TGK. In general, n000 did not coincide with n0 (GE2) or n00 (GE5). Therefore, it was convenient to differentiate the parameters n0, n00 and n000. In addition, the GE8 evidenced that n∞ depends on nobs, α and β, but not on n000 (n(t) = n∞ = n e for t → ∞). The parameters α and β in terms of u2, U1, θ, df, Df and nobs were given by

Eq (9) resulted from assuming that the value of n in the steady state was n = n e = (u2/U1)1/( and Eqs (7) and (8) were taken into account.

Simulations

Simulation of Eq (9)

Eq (9) coincided with Eq (6) for nobs = 1. The simulation of α (in days-1) versus df was shown for Df = 5 and four values for u2 (1, 10, 50 and 100 days-1) and nobs (1, 5, 10 and 20 cells). For this, values of df were varied from 0 to 5 with a step of 0.5, taking into account that df < Df. The simulation of β (in days-1) against df/Df was shown for four values of u2 (1, 10, 50 and 100 days-1) and the values of df/Df were ranged from 0 to 5 with a step of 0.5.

Simulations of GE2, GE5 and GE8

GE5 was used as reference because GE5 and GE8 were reported for the first time in the literature. The simulations of GE2, GE5 and GE8 were shown in a graph of n(t). Simulation of GE2 was made for different values of n0 (1x103, 1x104, 1x105 and 1x106 cells). Additionally, GE8 was simulated for three different situations: 1) nobs = 1 cell (GE8 and GE5 coincided) and different values of n00 (5, 10, 15, 20 and 25 cells); 2) nobs = 1x104 cells and different values of n000 (1x104, 5x104, 1x105 and 2x105 cells); and 3) n000 = 1x105 cells and different values of nobs (5x103, 1x104, 5x104 and 1x105 cells). In all these simulations, α = 1.0 days-1 and β = 0.3 days-1.

Experimental groups

In this study, V(t) was used by three reasons: 1) V(t) is related to n(t) and can be used interchangeably; 2) V(t) is less cumbersome to estimate than n(t) and it is frequently used in preclinical [9-11] and clinical [10] studies; and 3) the graphs of V(t) and n(t) shown sigmoidal shapes. Consequently, n(t) in GE1, GE2, GE5 and GE8 was replaced by V(t); n0 in GE2 by V0; n00 in GE5 by V00; n000 and nobs in GE8 by V000 and Vobs, respectively. In addition, nmed was replaced by Vmed and n∞ by V∞. The parameter V∞ was the tumor volume when t → ∞.

Two experimental groups were formed, each consisting of 10 male BALB/c/Cenp mice. The first group corresponded to the Ehrlich tumor, denominated G1, while the second group to the fibrosarcoma Sa-37 tumor, denominated G2. Experimental data of V(t) for Ehrlich and fibrosarcoma Sa-37 tumors were reported in [9], corresponding to their control groups.

Interpolation of experimental data

The Hermite interpolation method [13] was used to interpolate volume data of each individual tumor, in G1 and G2.

Estimation of values of α, β, df, Df and u2 from experimental data

Values of α and β (GE1, GE2, GE5 and GE8) and Vobs (GE8) were obtained from the individual fitting of each tumor volume (Ehrlich and fibrosarcoma Sa-37). The value of Vobs estimated directly with GE8 was named Vobs(α,β). The value V0 = V00 = V000 = 0.5 cm3 was the tumor volume chosen to describe TGK. This volume value was reached 15 days after 2x106 cells for the Ehrlich tumor and 5x105 cells for the fibrosarcoma tumor Sa-37 were inoculated in the BALB/c/Cenp mouse (see details in [9]).

Three equations in terms of df, Df and u2 resulted when Eq (6) was substituted in GE1, GE2 and GE5. The values of these three parameters were determined when each of these equations was used to fit experimental data. Besides, Eq (12) was substituted in GE8 and resulted an equation in terms of df, Df, u2 and Vobs, from which their values were estimated from fitting experimental data. Once known the values of df, Df, u2 and Vobs, they were substituted in their respective Eqs (6) and (9) to calculate their corresponding values of α and β. Values of α, β and Vobs obtained by this way were denominated αc, βc and Vobs(u2,df,Df), respectively, to distinguish these values from those that were directly obtained from fitting of the experimental data with GE1, GE2, GE5 and GE8.

The estimation errors for α, β, df, Df, u2, Vobs and Vobs(u2,df,Df) were denominated eα, eβ, edf, eDf, eu2, eVobs and eVobs(u2,df,Df), respectively. The estimation error for each parameter was reported for each individual tumor of Ehrlich and fibrosarcoma Sa-37.

The difference between α and αc, named Δα (Δα = α—αc), was calculated for each equation (GE1, GE2, GE5 and GE8) and experimental group (G1 and G2). In addition, it were computed differences between β and βc, denominated Δβ (Δβ = β—βc), and Vobs(u2,df,Df) and Vobs(α,β), denominated ΔVobs (ΔVobs = Vobs(α,β)—Vobs(u2,df,Df)).

Criteria for model assessment

Five quality-of-fit criteria were used for fitting of experimental data with GE1, GE2, GE5 and GE8: the sum of squares of errors, SSE (Eq (10)); standard error of the estimate, SE (Eq (11)); adjusted goodness-of-fit coefficient of multiple determination, (Eq (12)), that depended on goodness-of-fit coefficient r2 (Eq (14)); predicted residual error sum of squares, PRESS (Eq (14)); and multiple predicted residual sum error of squares, MPRESS (Eq (15)) [1, 3, 14], given by where was the j-th measured tumor volume at discrete time t, j = 1, 2, …, n1; was the j-th estimated tumor volume by GE1, GE2, GE5 and GE8; n1 the number of experimental points (n1 = 10) and k the number of parameters (k = 2 for GE1, GE2 and GE5, and k = 3 for GE8). The fitting was considered to be satisfactory when . Higher meant a better fit. was the estimated value of when GE1/GE2/GE5/GE8 was obtained without the j-th observation. MPRESS removed the last n1−m measurements. Each equation (GE1, GE2, GE5 and GE8) was fitted to the first m measured experimental points (m = 3, 4 or 5) and then from calculated model parameters the error between tumor volume estimated and measured values in the remaining n1−m points was calculated. Least Sum of Squares of Errors was obtained when SSE was minimized in the Marquardt-Levenberg optimization algorithm.

The Root Means Square Error, RMSE (Eq (16)) and the maximum distance, Dmax (Eq (17)) were also calculated following the methodology suggested in [1, 3, 14], given by where M was the number of interpolated data of tumor kinetics (graph of V(t)). Fi was the i-th tumor volume of the experimental data, which was chosen as reference. Gi was the i-th tumor volume calculated with GE1, GE2, GE5 and GE8.

Each fit with the GE1/GE2/GE5/GE8 was performed for each animal growth curve. A computer program was implemented in the Matlab® software (version R2012b 64-bit, Institute for Research in Mathematics and Applications, University of Zaragoza, Spain) to calculate the tumor volume. In addition, the mean ± standard error of each parameter of the equation (α, β, Vobs(α,β), u2, df, Df, Vobs(u2,df,Df), αc, βc), fit criterion (SE, PRESS, MPRESS, , RMSE and Dmax) and estimation error (eα, eβ, edf, eDf, eu2, eVobs and eVobs(u2,df,Df)) were calculated from their individual values, in each experimental group, following the methodology reported in [1, 3]. These calculations were performed on a PC with an Intel(R) core processor (TM) i7-3770 at 3.40 GHz with a Windows 10 operating system. All calculations took approximately 10 min, for each equation.

Results

Simulation of Eq (6)

Fig 1 showed the simulations of β versus df/Df (Fig 1) and α versus df (Fig 1) for different values of u2. The positive values of α (in the interval 0 ≤ df < 1) and β (in the interval df/Df < 1) increased non-linearly with the increase of df and decreased linearly with the increase of df/Df, respectively. The negative values of α increased non-linearly with the increase in df (df > 1.5). The negative values of β decreased linearly with the increase of df/Df (df/Df > 1). These behaviors were noticeable for the greater value of u2. Additionally, the parameter α had a discontinuity in the interval 1 < df < 1.5 and β = 0 when df/Df = 1 for all values of u2.

Fig 1

Simulation of Eq (6).

For different values of u2 (1, 10, 50 and 100 days-1) it is plotted (A) Graph of α (in days-1) versus df and (B) Graph of β (in days-1) versus df/Df.

Simulation of Eq (9)

Results of the simulation of β versus df/Df in Eq (11) coincided with that shown in Fig 1 (see Eqs (6) and (9)). The simulation of α versus df for nobs = 1 (Fig 2) reproduced the same result as in Fig 1. However, values of α were more negative, in the interval 0 ≤ df < 1, when nobs increased, being noticeable for the higher value of u2 (Fig 2). In Fig 2, as in Fig 1, it was observed a discontinuity of α in the interval 1 < df < 1.5.

Fig 2

Simulation of Eq (9).

For different values of u2 (1, 10, 50 and 100 days-1) it is plotted the graph of α (in days-1) versus df for (A) nobs = 1 cell. (B) nobs = 5 cells. (C) nobs = 10 cells. (D) nobs = 20 cells.

Simulations of GE2, GE5 and GE8

Fig 3 showed the behavior of n(t) when GE2 (Fig 3, GE5 (Fig 3) and GE8 (Fig 3) were used. Fig 3 revealed that the highest value of n∞ and the fastest TGK occurred for the highest values of n0 and α. Fig 3 showed that TGK was faster with the increase of n00 and all TGK tended to the same value of n∞ for all value of n00, keeping constant values of α and β. In this case, TGK was faster when the value of n00 increased with respect to nobs (Fig 3), being noticeable when nobs increased with respect to 1 (Fig 3). It is important to note that n0 = n00 (Fig 3) and n0 = n000 (Fig 3).

Fig 3

Evolution of the number of cells with time.

Simulation of the number of cells at time t, in days, (n (t)) for α = 1.0 days-1 and β = 0.3 days-1. (A) Simulation of GE2 for different values of n0 (1x103, 1x104, 1x105 and 1x106 cells). (B) Simulation of GE8 for nobs = 1 cell (coincides with GE5) and different values of n00 = n0 (5, 10, 15, 20 and 25 cells). (C) Simulation of GE8 for nobs = 1x104 cells and different values of n000 = n0 (1x104, 5x104, 1x105 and 2x105 cells). (D) Simulation of GE8 for n000 = n0 = 1x105 cells and different values of nobs (5x103, 1x104, 5x104 and 1x105 cells).

The results of Fig 3 showed that TGK grows slower (when n < n000) and then faster (when n > n000) for the greater value of nobs; all TGK were cut at t = 0 (same value of n000), for all value of nobs; and the value of n∞ depended on nobs and not n000 for each TGK. The results shown in Fig 3 were noticeable when the value of α increased with respect to that of β (results not shown).

Fitting of experimental data with GE1, GE2, GE5 and GE8 and estimation of values of α, β, df, Df and u2

The mean ± standard deviation of each parameter of the equation, fit criterion and estimation error were shown in Tables 1 and 2 of each equation (GE1, GE2, GE5 and GE8) used to fit experimental data of the Ehrlich and fibrosarcoma Sa-37 tumors, respectively. Tables 1 and 2 shown for these two tumor histological varieties: 0 < df < 1; 1 < Df < 2; 0 < u2 < 1; the highest values of α, u2 and the lowest values of df and Df for GE8; the lowest SE values for GE5 and GE8; the lowest values of PRESS, MPRESS, RMSE and Dmax; the highest values of r2 and for GE2, GE5 and GE8; and values of the parameter α differed when GE1, GE2, GE5 and GE8 were used. Nevertheless, the parameter β was the same when GE2, GE5 and GE8 were used, but not for GE1.

Table 2

Parameters of the models for the fibrosarcoma Sa-37 tumor.

Parameters	Different formulations of Gompertz equations
	GE1	GE2	GE5	GE8
α (days-1)	0.188±0.016	0.491±0.034	0.316±0.018	0.833±0.132
β (days-1)	0.127±0.017	0.252±0.018	0.252±0.018	0.252±0.018
Vobs(α,β) (cm3)	-	-	-	0.148±0.088
u2 (days-1)	0.274±0.093	0.530±0.152	0.471±0.132	0.576±0.070
df	0.759±0.074	0.822±0.070	0.746±0.058	0.688±0.042
Df	1.704±0.672	1.810±0.612	1.837±0.613	1.256±0.191
Vobs(u2, df,Df) (cm3)	-	-	-	0.142±0.029
αc (days-1)	0.197±0.020	0.494±0.029	0.322±0.011	0.814±0.082
βc (days-1)	0.152±0.018	0.290±0.020	0.280±0.017	0.252±0.018
SE	0.162±0.008	0.082±0.038	0.083±0.038	0.083±0.038
PRESS	0.761±0.227	0.063±0.059	0.063±0.059	0.064±0.060
MPRESS	0.623±0.203	0.063±0.059	0.064±0.059	0.064±0.060
r2	0.995±0.004	0.998±0.001	0.998±0.001	0.999±0.001
ra2	0.996±0.004	0.998±0.001	0.998±0.001	0.999±0.001
RMSE (cm3)	0.161±0.008	0.082±0.038	0.082±0.038	0.082±0.038
Dmax (cm3)	0.499±0.013	0.206±0.109	0.206±0.100	0.207±0.110
eα	0.025±0.011	0.046±0.022	0.061±0.012	0.079±0.035
eβ	0.034±0.009	0.053±0.013	0.057±0.029	0.055±0.023
eVobs(α,β)	-	-	-	0.027±0.007
eu2	0.031±0.003	0.035±0.013	0.039±0.010	0.061±0.015
edf	0.065±0.012	0.069±0.014	0.067±0.016	0.071±0.025
eDf	0.235±0.086	0.336±0.045	0.679±0.119	0.125±0.031
eVobs(u2,df,Df)	-	-	-	0.041±0.017

Means ± standard deviation of parameters of the fibrosarcoma Sa-37 tumor and criteria for model assessment obtained for different formulations of Gompertz equations.

For the Ehrlich tumor, Δα = 0.003, 0.005, 0.001 and 0.005 days-1 for GE1, GE2, GE5 and GE8, respectively. The variable Δβ = 0.012, 0.026, 0.014 and 0.000 days-1 for these respective equations and ΔVobs = 0.007 cm3. For the tumor fibrosarcoma Sa-37, Δα = 0.009, 0.003, 0.006 and 0.019 days-1 for GE1, GE2, GE5 and GE8, respectively. The variable Δβ = 0.025, 0.038, 0.028 and 0.000 days-1 for these respective equations and ΔVobs = 0.006 cm3.

Discussion

This study shows that GE2, GE5 and GE8 can be used interchangeably to describe experimental data of Ehrlich and fibrosarcoma Sa-37 tumors, taking into account their higher values of r2 and , and lower values of each parameter of the equation, fit criterion, estimation error, Δα, Δβ and ΔVobs (ΔVobs is only calculated for GE8).

The theoretical and experimental results of this work confirm different findings reported previously in the literature, such as: 1) the fractal origin of GE1, GE2, GE5 and GE8, as reported in [4, 15]; 2) the fractal property of tumors once reached nmed/Vmed, a matter that agrees with [16, 17]; 3) the role of the fractal dimension for the understanding of TGK, as suggested by Sokolov [18] and Breki et al. [19]; and 4) 1 < Df < 2, in agreement with [4, 20, 21] and the preferential growth along the largest diameter of the tumor, despite its ellipsoidal geometry [1, 3, 9, 11]. This fourth finding is in contradiction with 2 < Df < 3 reported by Breki et al. [19] in patients with metastatic melanoma; 5) The condition 0 < u2 < 1 for both types of tumors is consistent with the Steel equation [12]. If u2 = 0, then the tumor growth fraction must be high so that its mean doubling time (TD) is short, in contrast to [10, 12]. If u2 = 1 day-1 (all cancer cells are in apoptosis), TD → ∞ and α = 0 (tumor self-destruction), in contrast to the failure of the apoptosis mechanism in malignant tumors (because of the gene p-53 is repressed) and the existence of other cell loss mechanisms (metastasis, necrosis and exfoliation) [10, 11, 22]. The increase in u2 brings about a decrease in TD and therefore a higher value of α (Figs 1, 2).

Other novel findings have been revealed in this investigation that may be of interest for understanding of TGK, such as: 1) TGK sigmoidal form and n∞/V∞ do not depend on n0 and if on α, β and nobs/Vobs, when a given tumor histological variety grows in a certain type of syngeneic host to it. In this way, the action form of parameter n0/V0 (form or location) is eliminated in GE2, as reported in [2]. 2) The GE8 states that n0 in the GE2 is not a constant parameter but depends non-linearly with nobs/Vobs, n000/nobs (V000/Vobs), β and t. 3) The growth of a malignant tumor occurs for 0 < df < 1 and not when df = 0 (α = 0: the tumor does not form), 1 < df < 1.5 (discontinuity of α due to forbidden conformations or very unlikely tumor) and df > 1.5 (α < 0: the tumor self-destructs), in contrast to the values of df (1 < df < 2) reported in [4, 14, 23]. The forbidden conformations of the tumor can be explained by its stereochemistry due to the steric collides between all its elements and the tumor-host interaction. 4) The increase of α with the increase of df, at 0 < df < 1, confirms that the growth efficiency of a malignant tumor increases with its df, in agreement with [17, 24]. 5) Eq (11) states that this increase of α with df occurs if nobs satisfies strictly the condition ; otherwise, α < 0 for all β positive (Fig 2B, 2). The case α < 0 means that the tumor self-destructs, in contrast to the experiment.

The established condition for nobs suggests that: 1) nobs/Vobs depends on df and the ratio u2/β; 2) the fractal property of a malignant tumor also happens before or long before its detection (nmed/Vmed), as reported in [1, 25]; 3) the ratio u2/β may be an indirect indicator of the apoptosis-angiogenesis relationship reported in [26, 27]; 4) endogenous anti-angiogenic factors or inhibitors of angiogenesis (endostatin, angiostatin, among others) are present in the tumor before or long before reaching nmed/Vmed; 5) the term e− (see GE8 and the established condition for nobs) and the decrease of the parameter β with the increase of df/Df corroborate the essential role of angiogenesis process and the displacement of the balance between endogenous anti-angiogenic factors and endogenous pro-angiogenic factors towards these latter, when the tumor volume grows at time t, consistent with [10, 17, 22, 28, 29].

From the mathematical point of view, the condition 0 < df < 1 may suggest that the contours of Ehrlich and fibrosarcoma Sa-37 malignant tumors have zero area and/or they are totally disconnected. The first assumption confirms that these two types of tumors can be delimited from their surrounding healthy tissue, as in [9, 11]. The second hypothesis is based on proposition 2.5 [30]: “A set with dim F<1 is totally disconnected”. In this proposition, F is any set and dimH is the fractal dimension Hausdorff. It is important to note that, although the tumor boundary is wide, df < 1 if its only fractality is given by a totally disconnected line contained in that wide band.

From the biophysical point of view, the tumor contour totally disconnected can indicate the existence in it of pores/channels formed randomly of different sizes and shapes, changing in the time. This porous contour of a tumor may be related to the angiogenesis process (neo-formation of blood vessels), the formation of spicules by fragmentation of the contour into simple forms of molds (for example, triangles), roundness, irregular edge, anisotropy, roughness and compactness, findings reported in [1, 3, 10, 22, 31–34]. We believe that the tumor angiogenesis process can be regulated by the amount of pores/channels existing in its contour to interconnect with the surrounding healthy tissue. This hypothesis can corroborate that the angiogenesis of a malignant tumor is an emergency and regulated by the structural and conformational dynamic transformations that occur during TGK, as reported in [1]. On the contrary, if these pores/channels do not exist, the tumor would behave as an isolated system and would self-destruct, in contrast to the experiment.

Fig 3 deserves a careful interpretation, taking into account experimental results reported in the preclinical [1, 3, 9, 11, 14] and clinical [10] studies. The result of Fig 3 corresponds with the selection of different values of n0/V0 in the same TGK for different instants t0. For this case, in the experiment is guaranteed fixed co, cell viability, the tumor histological variety and the type of syngeneic host to it. The higher value of n0/V0 in the same TGK means a larger tumor size, which is reached at a higher t0.

Results of Fig 3 are associated to the same tumor histological variety that grows in several types of syngeneic hosts to it. For this case, co and cell viability fixed are guaranteed, taking into account the role of the immune system in the delay of TGK, depending on its immunocompetence degree [10, 11, 22, 35]. As a result, tumors reach different values of n00/V00 o n000/V000 at the same time t0. The higher value of n00/V00 (n0/V0 in Fig 3) or n000/V000 (n0/V0 in Fig 3) corresponds to the lower immunocompetence degree of the host (e.g., an immunosuppressed host).

Results of Fig 3 refer to two possible situations: 1) different tumor histological varieties that grow in the same type of syngeneic host to them. For this case, co is different so that each tumor histological variety reaches the same value of n000/V000 at the same time t0. 2) A given tumor histological variety that grows in different types of syngeneic hosts to it. For this case, co is the same for each tumor histological variety. For these two cases, nobs/Vobs for each tumor histological variety is reached in a different tobs, in accordance with the experiment [9, 11]. These two situations become noticeable when β approaches α (results not shown). Furthermore, this figure reveals that for the highest value of nobs/Vobs (reached in a greater tobs) TGK is slower for n(t) < n000 (V(t) < V000) and then faster for n(t) > n000 (V(t) > V000). By contrast, the tumor that has the lowest nobs/Vobs is the fastest growing for n(t) < n000 (V(t) < V000) and then its TGK is slowest for n(t) > n000 (V(t) > V000).

The advantages of GE8 over the various formulations of GE [2, 3], the Hahnfeldt model [36-38] and mKJMA equation [1], used to describe undisturbed TGK, are: 1) inclusion of two parameters (nobs/Vobs y n000/V000) that are measured and estimated from experimental data. 2) TGK and n∞/V∞ can be known a priori if nobs/Vobs (starting point of TGK), reached at tobs, is estimated for each type of tumor that grows in a syngeneic host to it, as reported in [1, 3, 11].

The relation of the tumor growth with df and Df is previously obtained by using a mesoscopic formalism and fractal dimension [39]. Besides, Izquierdo-Kurlich [39] report the differences between df and Df and propose a relation between df and the dynamic quotient on the interface, named kc, (see Eq (48)). This relationship differs from that reported in [4] (see Eq (3)), which is used to obtain Eq (8). If the relation published in [39] is taken into account in this study, Eq (8) is also obtained, except a small change in α numerator (1/2 instead of 1). As a result, 0.75 and 1 are the discontinuities of α, instead of 1 and 1.5, respectively. Nevertheless, these change do not affect significantly the results of this manuscript and confirm that tumors exits for 0 < df < 1. It can be verified that df for Ehrlich and fibrosarcoma Sa-37 tumors are less than 0.75 and 1 when Eq (48) in [39] and Eq (3) in [4] are used.

In this study, the tumor growth in the time results of the complex interactions that happen in the tumor and between it and the surrounding healthy tissue, as in [3,14]. Nevertheless, in it does not explicitly discuss the interactions among the individuals neither the cooperative capacity of they in a population to explain its growth behavior, as in [25, 5–8]. These works confirm the fractal property of the tumors, as in this study. Therefore, an additional study may include these interactions for Eq (8).

Further studies can be carried out to validate GE8 in TGK of different tumor histological varieties that grow in both immune-competent and immune-deficient organisms. This will allow us to know how Df, df, u2, Vobs(α,β) and Vobs(u2,df,Df) change when using different types of tumors and degrees of immune-competence of several organisms, as well as confirming the relationship of these five parameters with the aggressiveness [1], angiogenesis [17], coherence [15, 16], anisotropy, heterogeneity, hardness, changes in the mechanical-elastic-electrical properties of a tumor, among others findings [1].

Conclusions

GE8 describes well the growth kinetics of the Ehrlich and fibrosarcoma Sa-37 tumors and includes two parameters that are directly estimated from the experiment that confirm the fractal property of the tumors and the fractal origin of different Gompertz formulations.

Appendix A

In [4] it is assumed that the growth ratio of the number n(t) of tumor cells obeys to the differential equation where m represents the number of tumor cells at the boundary of the tumor, u1 is the constant of the velocity of the mitosis and u2 is the constant of the velocity of apoptosis.

Assuming that the boundary has a fractal structure with dimension d, then , r being the average radius of the tumor. On the other side, n depends on the morphology of the tumor, described by the fractal dimension D, and . The morphological constants k1 and k2 are related to the magnification of the image [4].

Substituting these values of m and n and eliminating r, Eq (1) can be written as a Bertalanffy-Richards equation. where n = (u2/U1)1/(1− is the value of n at the steady state, the dimensionless morphological parameter θ is defined by θ = d/D and U1 is given by .

Taking into account that the above equation is approximated in [4] by the Gompertz equation

This approximation is valid when θ→1 or n→n.

In [36] it is justified that the quotient U1/u2 can be expressed as a function of d and in [4] it is shown that the solution of the differential system (2) can be expressed as a Gompertz equation (Eq (1) in this paper) with the intrinsic growth rate of the undisturbed tumor, named α (α > 0), and the deceleration factor, named β (β > 0), related to the tumor fractal dimensions by

Supporting information.

(TXT)

Click here for additional data file.