Gompertz models with periodical treatment and applications to prostate cancer

1.   Introduction

Prostate cancer (PCa) is the second-most frequently diagnosed cancer and the second-most frequent cause of cancer death in men ^[1,2]. It was estimated that approximately 1.2 million new cases of PCa are diagnosed and 350,000 men die of the disease each year ^[3]. Globally, the incidence of PCa is increasing, mainly among developing countries. The American Cancer Society estimates that in 2023, there were about 288,300 new cases of PCa and 34,700 deaths in the United States ^[4].

Treatment of PCa varies based the stage of the cancer. Those with localized disease at low risk for spread are monitored regularly by active surveillance. Those at higher risk may receive treatment (prostatectomy or radiation therapy) to eliminate the tumor. Those with metastatic disease are treated with chemotherapy or other agents to alleviate the symptoms of metastatic tumors ^[5]. Despite having a regression initially in most cases, the PCa usually evolves into castration-resistant prostate cancer (CRPC), which has a lethal result. In this scenario, docetaxel is the preferred chemotherapeutic agent available during the treatment ^[1,6]. However, the patient's survival rate still depends on multiple factors, among those is the acquired drug resistance, which is a clinical issue, especially for docetaxel. In this regard, combination therapy proves to be more efficient and is becoming quite relevant among PCa therapies ^[1,2,6]. Specifically, the combination treatment performed with docetaxel and capsaicin (natural nutrition compound extracted from hot chili peppers) blocks tumor growth in vivo efficiently as demonstrated in ^[7].

Mathematical modeling is useful for quantitative descriptions of plenty of physiopathological phenomena by generating valuable clinical foresight. Due to a lack of knowledge about the complexity of the PCa, which is also related to several complex biological mechanisms that take place during different types of treatment, it becomes essential to include mathematical models to better explain the dynamics of the interactions of all variables in these physiopathological phenomena. In the last two decades, many mathematical models for prostate cancer have been developed, primarily because of the amount of experimental and clinical data collected due to technological achievements and the success of the earliest models. We refer to Phan et al. ^[8] for a review on various mathematical models of prostate cancer and clinical applications.

For many decades, tumor growth kinetics has been studied experimentally and clinically ^[9]. A relevant finding is that tumor growth is only partially exponential since the growth rate decelerates after a certain amount of time. The most widely accepted cancer model is the Gompertz model which has been used in numerous studies involving human data ^{[10,11,12,13,14]} and takes the form:

where $r, K$ are positive constants respectively representing the growth in proportion to cell population size and the maximal volume (carrying capacity), $V_{0}$ is the initial tumor volume, and $t \geq 0$ is the time, as first described in ^[10].

Note that treatment is usually periodic ^[1,6]. In this paper, we propose a Gompertz model with periodic treatment to mathematically represent cancer tumor growth under periodic therapy, taking into account anticancer drug's saturation. First, we establish the existence, uniqueness, and stability of a periodic solution of the model. Second, we apply our mathematical model to calibrate the data reported in ^[7] on the treatment of PCa. Our simulations confirm the conclusion of ^[7], that is, the combination treatment performed with the union between the two anticancer drugs, capsaicin and docetaxel, is the most efficient treatment, among the ones analyzed.

2.   Gompertz model with periodic treatment

Our goal is to understand temporal tumor volume behavior when a periodical treatment is provided. As shown in ^[10], the growth of a tumor can be described by a Gompertz curve. In other words, if $x(t)$ denotes the volume at time $t$, then the dynamics of $x(t)$ are governed by the following differential equation

where $r$ is the growth in proportion to cell population size and $K$ is the maximal volume (carrying capacity). As mentioned previously, the Gompertz model (2.1) is a commonly used empirical description for tumor growth. Other models, such as the logistic and exponential equations are also utilized ^[15,16]. Despite the widespread use of these models in explaining experimental data, a deeper understanding of their success is still required. A comprehensive model that derives these three models from specific biophysical assumptions concerning various types of cellular interactions linking microscale variables related to cell proliferation to the macroscale variable population can be found in ^[17].

To introduce the effect produced by the treatment, we propose the following model:

where $h(t)$ is a $\omega$-periodic function that mimics the periodic dosage, representing the proportion of tumor killed at time $t$. For this, we assume that $0 \leq h(t) \leq 1$. Hence, the term $h(t)x(t)$ in (2.2) can be interpreted as the induced volume decrease at $t$.

Following ^[18] and ^[19], we assume that the function $h(t)$ takes the following form:

with $F$ and $c$ being positive constants. We will see that, although model (2.2) does not reflect the actual behavior of the volume when compared to the data, it gives us relevant results about the global stability of a more general model considering saturation and resistance.

3.   Periodic treatment with saturation

Now consider the case where treatment effectiveness changes over time. The treatment process is called saturation if the drug gets more efficient until it reaches a plateau. To include this phenomenon in the dynamics of $x(t)$, we consider the following modification of (2.2):

where $s(t)$ is a function describing the behavior of saturation.

The Hill equation is widely used to model the pharmacologic effect of drugs ^[20,21] and is often expressed in the following way ^[22]:

where $A$ is the drug concentration that produces a $50\%$ maximal response, $C$ is the total drug concentration at time $t$ and $n$ is called the Hill constant.

The total drug concentration $C(t)$ can be written as ^[23]:

where $P$ represents the highest level of drug in the bloodstream during the distribution phase, while $q$ indicates the entry rate of drug into the tumor from the blood. It is important to note that $q$ is dependent on $\omega$. By substituting (3.3) into (3.2), we obtain

Taking $D = (A/P)^n$ and $\lambda = qn$, the Hill Eq (3.2) turns into the logistic function

Thus, as a prototype for the saturation function, we propose

Notice that $s(0) = 1/(D+1) < 1$ corresponds to the initial fraction of the real effectiveness of the drug at the beginning of the treatment. Moreover, $s(t) \rightarrow 1$ as $t \rightarrow \infty$, meaning that the drug approaches its maximum effectiveness as $t$ increases.

4.   Periodic solutions and global stability

In this section, we will prove that (2.2) has only one periodic solution, which is globally stable.

Theorem 4.1. The Gompertz model (2.2) with periodic treatment $h(t)$ given by (2.3) has a unique periodic solution which is globally stable.

By applying the change of variables $y = \ln x$, we know that (2.2) is topologically equivalent to

for $h$ given in (2.3). Notice that there is a one-to-one correspondence between periodic solutions of (4.1) and (2.2). Indeed, if $\Phi(t)$ is a periodic solution of (4.1), then $\exp(\Phi(t))$ is a periodic solution of the original Eq (2.2). Similarly, if $\Psi(t)$ is a periodic solution of (2.2), then $\ln(\Psi(t))$ is a periodic solution of (4.1). Therefore, it suffices to prove the existence and uniqueness of $\Phi(t)$ for Eq (4.1).

4.1.   Existence

The solutions of (4.1) are of the form:

For all $y_0$, the solutions (4.2) are bounded. Indeed, using the fact that $h$ is a product of bounded functions, and hence bounded, we obtain that

Thus, By Theorem $4.11$ in ^[24], Eq (4.1) has an $\omega$-periodic solution, namely, $\Phi(t)$. It is noteworthy that, although the theorem is originally stated for $1$-periodic functions, the result holds true for $\omega$-periodic functions in general. This can be verified by rescaling the time, which is explained in Section $4.2$ of ^[24].

4.2.   Uniqueness

Suppose that there exist two different $1$-periodic solutions of (4.1), $\omega(t)$ and $\Phi(t)$. Since their graphs do not intersect, we may consider, without loss of generality, that $\omega(t) > \Phi(t), \ \forall t \geq 0.$ Now, notice that the curve $v(t) : = \omega(t) - \Phi(t)$ is a $1$-periodic curve into $\mathbf{R}_{+}$. Moreover,

Thus, for all $t \geq 0$ we have

Integrating the above expression from 0 to 1, we get that

since $r > 0$, contradicting the periodicity of $v$.

4.3.   Stability

We can compute the Poincaré map $P(y) : = \varphi(1, 0, y)$, given by

Let $y \in \{ \Phi(t); t\geq 0\}$. Thus,

Hence, by Lemma 4.21 and Theorem 4.22 of ^[24], $\Phi(t)$ is an asymptotically stable periodic solution.

Consider $y_{0} : = \Phi(0)$, then $y_{0}$ is a fixed point of the Poincaré map, i.e.,

Thus, by (4.3), we get

Therefore, we are able to have the explicit formula for the 1-periodic solution

due to (4.2).

Now we observe that not only the $1$-periodic solution $\Phi(t)$ is asymptotically stable, but it is actually globally stable. Indeed, by (4.2), given another initial value $x_{0} > 0$, we have that

and, hence,

proving the global stability of $\Phi(t)$.

A direct implication of this fact is that (2.2) has a unique periodic solution $\Psi(t)$, which is $1$-periodic, globally stable, and given by

4.4.   Model (3.1) with saturated treatment

Notice that $s(t) = \frac{1}{De^{-\lambda t}+1}$ given in (3.4) is time-dependent but not periodic, so model (3.1) with saturated treatment is no longer time periodic. By taking $y = \ln(x)$ we have

The solutions of (4.4) are of the form:

Using the expressions of $s(t)$ and $h(t),$ we further have that

Note that we have the following estimates:

By taking the limit as $t\to \infty$ in Eq (4.6), we obtain the upper and lower bounds of the limit of the trajectories.

Proposition 4.2. For all $y_0\in \mathbb{R}$, we have

5.   Numerical simulations and data fitting

In ^[7], the authors studied the efficacy of the combination of docetaxel and capsaicin on prostate cancer cell proliferation. The study was conducted over 15 days, and the analyzed treatments consisted of daily i.p. injections of Vehicle dimethyl sulfoxide (DMSO), $2$ mg/kg capsaicin (CAP), $10$ mg/kg docetaxel (DTX), or $2$ mg/kg capsaicin $+ 10$ mg/kg docetaxel (CAP + DTX).

Using the data available in ^[7] for LNCaP (Lymph Node Carcinoma of the Prostate) cells, we will fit both models:

and

Here (5.1) and (5.2) are particular cases of (2.2) and (3.1), respectively. Since the data is from a short-term experiment, we did not include the resistance term in (5.2). Moreover, $\omega = 1 \ day^{-1}$ once the drugs were administrated daily.

All parameters were calculated using the least squares function in Python from scypy.optimize package.

5.1.   Model (2.1): Standard Gompertz model

By utilizing the control group, we successfully fitted the model (2.1) and subsequently computed the parameters $r$ and $K$ (Table 5.1). Since these constants represent intrinsic characteristics of the tumor, they must be the same as the ones appearing in (5.2) and (5.1). Unless explicitly specified otherwise, all plots share the common initial condition $x_0 = 83.61\mathrm{mm}^3$, corresponding to the average initial volume of tumors analyzed in ^[7].

5.2.   Model (5.1): Gompertz with periodical treatment

Having fixed the values obtained for $r$ and $K$ in Table 5.1, we proceed to apply the same method as before, but this time in order to estimate the values of the parameters $F$ and $c$ that correspond to the treatment component of the model (5.1). The resulting parameter values are presented in Table 5.2 and the corresponding fit for each treatment protocol is depicted in Figure 5.2.

5.3.   Model (5.2): Gompertz with saturation

Repeating the process for model (5.1), wherein we still fix the values for r and K in Table 5.1, but now we estimate the values of the parameters $F$, $c$, $D$, and $\lambda$ appearing in the treatment component with saturation of the model (5.2). The resulting parameter values are presented in Table 5.3, and the corresponding fit for each treatment protocol is depicted in Figure 5.3.

5.4.   Periodic solution and global stability

Having obtained all parameters, it is now feasible to graphically represent the globally stable periodic solution of model (5.1). This is achieved by employing the initial condition stipulated using the Poincaré map (4.3). To elucidate, we numerically computed the global stable periodic solution for the model with parameters corresponding to the CAP treatment (see Table 5.2). This procedure is applicable analogously to the remaining two configurations. In Figure 5.4, we depict the periodic solution starting at $x_0$ = $9.98$ $\mbox{cm}^3$, together with proximate solutions, thereby illustrating its inherently attractive behavior.

5.5.   Model (5.1) vs Model (5.2)

Both models (5.2) and (5.1) have effectively captured the dynamics of the tumor volume during capsaicin treatment, as demonstrated in Figures 5.3(a) and 5.2(a). Upon comparing the parameters obtained for both models, it is clear that the values assigned to $F$ and $c$ exhibit a remarkable similarity under the CAP treatment (see Tables 5.2 and 5.3). Furthermore, the parameter $D$ in (5.2) is close to zero, which signifies that capsaicin treatment does not induce substantial saturation effects in this context. This explains why both models agree with each other, producing similar results.

Nevertheless, in instances involving the use of docetaxel, model (5.2) distinguishes itself by precisely aligning with the data (refer to Figure 5.3(b)), affirming the occurrence of the saturation process. This phenomenon could be attributed to factors such as dosage, administration methods, or the distinctive pharmacokinetic characteristics inherent to the administered drug ^[25].

5.6.   Long Term Trajectories

In Figure 5.5, we plotted the solutions in the interval from $t = 0$ to $t = 100$ with initial value $x_0 = 83.61$ $\mathrm{mm}^3$. Corroborating the conclusion in ^[7], we see that the docetaxel and capsaicin synergistic effect is the most efficient in reducing tumor growth when compared with the other two protocols (see Figure 5.5(c)).

In Figure 5.6, we extended the temporal range from 0 to 400, placing emphasis on the trajectories and their convergence toward the region delineated by Proposition 4.2.

6.   Discussion and Conclusions

In the present study, we have developed and validated a comprehensive framework that characterizes tumor size dynamics under periodic treatment. Our models provide an accurate representation of the therapeutic effect in the case of prostate cancer, enabling us to forecast the evolution of treatment.

Our study introduces two models, one being a more general version of the other, encompassing the saturation effect during treatment. By analyzing the dynamics, we demonstrated the global stability of both models and established the existence and uniqueness of a periodic solution for (5.1). Furthermore, we illustrated the efficient fit of both models to capsaicin treatment data for prostate cancer. Our findings indicated minimal saturation during capsaicin treatment, and model (5.1) aptly described dynamics in treatments where saturation was not observed.

Aligning with conclusions in ^[19], our plotted solutions revealed that a combined treatment of capsaicin and docetaxel is the most efficient among the analyzed treatments. Importantly, our model suggests that this treatment regimen is the only one that potentially leads to complete remission.

Beyond PCa, our models have the potential to be useful in the analysis of other cancers. Similarly, we anticipate the potential modification of the model to incorporate drug tolerance, a crucial consideration given the known role of genetic mutations, phenotypic plasticity, and stochastic cell-to-cell variability in driving drug resistance and tumor relapse ^[26]. While acknowledging the relevance of these phenomena, validating our model in this context necessitates data from more extensive experiments, representing an avenue for future analysis.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

We thank two anonymous reviewers for their helpful comments and suggestions.

Conflict of interest

The authors declare there is no conflict of interest.