Abstract
Background/Aim: Mathematical models have long been considered as important tools in cancer biology and therapy. Herein, we present an advanced non-linear mathematical model that can predict accurately the effect of an anticancer agent on the growth of a solid tumor. Materials and Methods: Advanced non-linear mathematical optimization techniques and human-to-mouse experimental data were used to develop a tumor growth inhibition (TGI) estimation model. Results: Using this mathematical model, we could accurately predict the tumor mass in a human-to-mouse pancreatic ductal adenocarcinoma (PDAC) xenograft under gemcitabine treatment up to five time periods (points) ahead of the last treatment. Conclusion: The ability of the identified TGI dynamic model to perform satisfactory short-term predictions of the tumor growth for up to five time periods ahead was investigated, evaluated and validated for the first time. Such a prediction model could not only assist the pre-clinical testing of putative anticancer agents, but also the early modification of a chemotherapy schedule towards increased efficacy.
- Pharmacokinetic (PK)–Pharmacodynamic (PD)
- tumor growth inhibition (TGI) mathematical model
- deep learning neural networks (DLNN)
- nonlinear optimization
- TGI model parameters estimation
- adaptive tumor growth short-term prediction
- xenografted mice (PDX)
- pancreatic ductal adenocarcinoma (PDAC) xenograft
Despite the major scientific advances in cancer biology and treatment during the last decades, the disease remains the second leading cause of death according to the World Health Organization (WHO) (1). In order to understand the response and detect the progression of a tumor in relatively short time without the huge costs of clinical and laboratory experiments, tumor growth mathematical modeling has been extensively developed and used during the last decades. Mathematical models have been considered as tools towards understanding cancer biology and aiding therapy since the 1950's (2). Classical mathematical models such as the logistic Gompertz and von Bertalanffy have been extensively tested and used due to their ability to predict the course of tumor progression (3). Indeed, several models, based on nonlinear mathematical equations and on tumor growth physiology and biology have been proposed as tools to predict the growth of a tumor (4-6). Multiscale tumor growth modeling, dynamic differential equations models, delay differential mathematical models, etc., (7-9), are capable of modeling more efficiently the multiscale and the spatio-temporal evolution of tumors. Although these models are capable of simulating and computing tumor cells' growth accurately, some of them incorporate the inputs not explicitly and in a very complex way, i.e. time delayed, or they do not incorporate model inputs at all. Thus, these models are not able to correlate explicitly the anticancer treatment dosing schedule, as input, to the observed tumor growth, as output. The later though is necessary as it consists the only approach that allows us to take into consideration quantitatively and for the total treatment time period the effects of the administered chemotherapy drugs in the tumor growth evolution. Thus, a significant problem generally encountered in all the existing mathematical models that deal only with the tumor growth is that this fact prohibits ultimately any possibility of their use for direct assessment of the anticancer effect of an agent with computational and numerical studies.
A simple dynamic mathematical model, which is based on a small number of biologically and drug relevant parameters, linking the growth dynamics of a malignant neoplasm with the administration of an anticancer drug and thus, capable of simulating tumor growth inhibition with accuracy during treatment, is the Tumor Growth Inhibition (TGI) mathematical model introduced by Simeoni et al. in 2004 (10, 11). This TGI dynamic mathematical model has been extensively used for ranking anticancer candidates for their efficacy and for designing in vivo experiments successfully. However, its ability to predict the future tumor growth based on experimental data from a previous time period has not been examined yet (12, 13). In this context the aim of the present study was to identify and investigate the predictive power of such an important input-output dynamic nonlinear mathematical model-computational tool, by first estimating accurately its unknown PK-PD principal parameters based on (i.e. using) the in vivo tumor growth under anticancer agent treatment experimental data for the total time period of the treatment or for any part of it.
Materials and Methods
Nonlinear dynamic mathematical model of the tumor growth inhibition (TGI). The pharmacokinetic-pharmacodynamic (PK-PD) state-space mathematical model (TGI model, in brief) first presented by Simeoni et al. in 2004 (10-12), consists a nonlinear first-order differential equations mathematical model of the form: Ż(t)=f(Z(t), x)+B (Z(t))u(t), with four state variables Zi(t), i=0,1,2,3, input u(t) ≜ c(t)(c(t): the concentration of the administrated anticancer agent in mg/L) and a vector of five biologically relevant principal parameters xT ≜ [k1,k2,λ0,λ1,w0].
Human-to-mouse xenograft experimental data. The experimental data from human-to-mouse cancer xenografts used in the present study were obtained from the experimental studies reported by Bilalis et al. (14). In short, AsPC-1 human pancreatic cancer cell line was xenografted subcutaneously at the rear flank of male mice. All mice received intraperitoneal (i.p.) injections of gemcitabine twice in a week interval (days 19 and 26 after the inoculation of the tumor cells). The dose of the chemotherapeutic drug was at 100 mg/kg/injection/mouse. Tumor volumes were converted to masses according to the formula:
[1]
where ϱ=1 mg/mm3 is assumed to be the tumor tissue density (10, 14). The pharmacokinetic (PK) data of gemcitabine were obtained and depicted from the almost identical experiments found in the literature, i.e. Veerman et al. (15). One compartment model was also used to describe the pharmacokinetics of the administrated anticancer drug gemcitabine.
The estimation problem of the initially unknown parameters of the TGI nonlinear mathematical model. In many mathematical modeling cases of physical processes, the mathematical model parameters are unknown, and they must be estimated using the available input-output experimental observations. In the case of the experimental data of Bilalis et al. (14), the tumor growth curve (i.e. the TGI model output) in combination with the delivered chronological gemcitabine dosing (i.e. TGI model input) are used in order to estimate the TGI dynamic model's principal parameters values from the experimental observations (the physical process). More specifically, an objective functiondepending on the TGI model's unknown parameters values and the time, created by the comparator procedure of Figure 1, must be minimized with respect to the vector of the unknown parameters xT=[k1,k2,λ0,λ1,w0] ≜ [x1,…, xn], n=5 in the present case. For any given set (i.e. vector x) of the unknown parameters' values, the nonlinear TGI mathematical model, when solved numerically for the same time period as that of the experiments, simulates the tumor growth for the same time period, i.e. simulates chronologically the experiment. Thus, using the experimental tumor growth curve, the function is evaluated by the eqs. [2], [3] below, where E(x, t) is defined as the difference between the observed (measured) tumor mass at each time point that a measurement takes place and the tumor mass estimated by the TGI model simulation. The vector of the unknown parameters' values that minimizes the cost function J(x, t) of the eq. [3] is the solution to the TGI mathematical model's unknown primary parameters k1, k2, w0, λ0 and λ1 estimation problem.
Comparator procedure between the physical experimental process of xenografted mice and the under fitting to the experimental data mathematical TGI model.


In that way an accurate dynamic first order differential equations set able to simulate (i.e. to calculate numerically by integrating them) accurately the specific experiment is obtained.
Advanced non-linear mathematical optimization techniques, based on recent deep learning neural networks (DLNN). The DLNN adaptive unsupervised nonlinear constraint optimization technique proposed by Sadollah et al. (16, 17) searches to achieve the global minimum value of an objective function J(x1, …, xn, t) subject to constraints of the form xr,min≤xr≤xr,max, r=1, …, n.

Each “individual” defined as a set containing a numerical value of each optimization variable (unknown parameter) is called “pattern solution” (e.g., in the Genetic Algorithms this vector is called “Chromosome”). In a n dimensional optimization problem, a pattern solution is a vector of 1×n, representing input data in the DLNN. This 1×n “array” is defined as the xT above. To start the optimization algorithm, a candidate of pattern solution matrix with size Npop×n is generated. Hence, matrix which is randomly generated in the range between lower and upper bounds of each unknown variable (i.e., lower and upper bounds are assumed to be defined by a decision maker), is given as follows (rows and column are the population size (Npop) and the dimension size (n), respectively):
[5]
Next, their objective function J(x1, …, x5, t)=J(x1, …, x5,Z(t), u(t)) is evaluated following the procedure described in Figure 1 and eqs. [2], [3] and the initial “target” solution xtarget (i.e. the vector supposed initially that corresponds to minimum J value) is determined. Similarly, to the population of pattern solutions, a weight matrix W is generated:
[6]
The weights are generated randomly following uniform distribution and must satisfy the constraints imposed by the eq. [7] below. The existence of such constrains is particularly important as it prevents the algorithm from getting “stuck” at a local minimum (16).

Based on the previous population of pattern solutions and inspired by the weight summation technique usually used in ANNs (artificial neural networks) a new updated X population is generated and the weighting matrix W is updated using the equations below:
[8]
[9]
[10]
where k is the optimization iteration index. To prevent the algorithm from an early convergence, a bias strategy (by an operator β) is introduced (16). A certain fraction of the pattern solutions and weights generated in the new population is modified, making the bias operator to act as a noise of the system and providing the ability to widen the search space and explore other places, which have not been visited by the population yet. The probability to modify and generate new pattern solutions using the bias operator is decreased as the iteration index k is increased. Instead, it is more possible for a solution to be transferred towards the best solution (i.e. new “target” solution), see eq. [11] below:
[11]
The new “target” solution and its weight are next determined, and the above update procedure is repeated until the maximum number of iterations kmax is reached. The parameters estimation procedure is then restarted from the converged solution of the previous trial, thus a “rough check” is performed to see if the solution obtained corresponds to the global minimum of the objective function. In Figure 2, the flowchart of the TGI model's unknown parameter estimation procedure is presented. Table I presents the lower and upper physically imposed bounds of the TGI mathematical model's primary parameters values in order to achieve the best fit to the experimental data.
Adaptive measurement and short term ahead TGI prediction update An innovative methodology is proposed herein in order to perform statistically accurate forecast and evaluate the ability to accurately predict the tumor evolution in the near future after each time point te based on the experimental data (observations) available until that time instant te. The proposed methodology can be understood through the simple flow diagram shown in Figure 3.
The DLNN technique described above is first applied to estimate the TGI nonlinear mathematical model's principal parameters (i.e. k1, k2, w0, λ0 and λ1 and based on a subset of the available experimental data [0, te], te≤tend, tend, being the last day of the xenograft experiment. Then, the fitted TGI model is integrated numerically (i.e. simulated) with inputs the anticancer agent dosing of the experiment delivered until the day te (model inputs) for the time periods [0, te+1], [0, te+2], [0, te+3], [0, te+4], [0, te+5] in order to predict (calculate) the future tumor growth inhibition that could be achieved at a horizon ofor days (time periods) ahead, respectively. In the next step, the real measured tumor mass value at the day te+1 is subtracted from its corresponding predicted as above value and the prediction error norm (i.e. the prediction error absolute value) and its root mean square are calculated for the one, two, three, four, five time periods ahead predictions, respectively. Then, the update of the TGI model's parameters estimation is performed: the TGI model's parameters are re-estimated for the “new” updated time period [0, te,new=te+1] using the “new” extended subset of experimental data available until the day te,new=te+1. The time instant te,new is assigned to the te “time index” and the procedure re-starts, repeated as above and terminated when te,new=tend. However, in order to achieve more accurate predictions with lower errors and less computational efforts at the same time, a moving (rolling) window of appropriate length was applied on the experimental data as it is shown in the flowchart of Figure 3. A moving window is a way to isolate subsets of a long string of time-dependent measurements, simply by taking at each time instant the last l time instants measurements (t=n•ΔT, n=0,1, … N and l << N). An optimal selection of the moving window's length l=5 was achieved by applying the Akaike's information criterion (AIC) (18).
Flow chart of the DLNN adaptive unsupervised nonlinear constraint optimization algorithm for the estimation of the unknown parameters of the TGI mathematical model.
By updating the TGI model parameters estimation using the 5 “input-output” measurements of the moving window as it is formed at each time instant, an efficient track of the varying slopes of the experimental data curves can be achieved. In Figure 4, an explanatory graphical example of the above adaptive tumor growth short term ahead prediction procedure to assess the efficacy of an anticancer agent is shown [see also (19)].
Value range of the TGI mathematical model principal parameters.
Results
The TGI non-linear mathematical model principal parameters estimation. The pharmacokinetics of gemcitabine in mice were described using one-compartmental model, thus in Figure 5 the in-silico plasma concentration of gemcitabine after two i.p. injections at 100 mg/kg (in a 7-day interval) is shown (15). Pharmacodynamics after two i.p. injections of 100 mg/kg on days 19 and 26 after the inoculation were experimentally registered and are given in the Table II.
Linear interpolation was used to estimate the unobserved tumor masses in the time instants that an observation (tumor measurement) was not performed. The tumor mass average observed/interpolated during the experiments, and the TGI model tumor masses obtained following its numerical integration during the same, as in the actual experiment, time period, are shown in Figure 6. As it is shown in Figure 6, the identified TGI model could describe with high accuracy the tumor growth experimental curves (see the prediction's statistical errors in Table III). Using the estimated TGI model's primary parameters values (Table IV), a threshold concentration for the tumor eradication of ~336 ng/ml (0.336 mg/l) was estimated, which means that any drug concentration over this value eventually may lead to a delayed tumor growth and hence can be considered effective (10). The tumor time efficacy index (TEI) value, i.e. the efficacy of a treatment measured using the achieved delay of the tumor growth, was estimated to be ~7.2 days which means that the tumor growth was inhibited in this experiment by gemcitabine for about 7 days (10).
Adaptive short term ahead prediction of Tumor Growth Inhibition by administrating gemcitabine. Using the DLNN algorithm described above, the TGI model was first identified using the experimental input-output data until the 25th day (te=25, initially) (see Point 1 in Figure 4), then it was integrated for one, two, three, four, five time periods ahead in order to predict the tumor's growth at the 26th, 27th, 28th, 29th, and 30th day respectively (see Point 2 in Figure 4). On the 26th day, a dose of 100 mg/kg gemcitabine was administrated in mice. An optimal length of the moving window was applied based on the Akaike's information l=5 criterion (i.e. the last l time instants measurements). Thus, the TGI model parameters were re-estimated using the TGI data from the 22nd day to the 26th day of the experiment, i.e. te,new=te+1=26. Then, the prediction error was calculated by subtracting the predicted value from the corresponding real value, being known by the laboratory experiments. Starting from the new te=te,new=26, the tumor growth was predicted again for one, two, three, four, five time periods (days) ahead, i.e. for the days t=27, 28, 29, 30 and 31. Subsequently, the model's parameters were re-estimated using the TGI data from the 23rd to the 27th day of the experiment, i.e. te,new=te+1=27 (moving the window - see Point 3 in Figure 4) and the prediction error was calculated. The whole procedure of the measurement update, of the re-estimation of the TGI model's parameters following each measurement update (based on the last l=5 time instants measurements) and of the adaptive tumor growth prediction is continued step by step (see the flowchart of Figure 3 and the predictions in Figure 4). Interestingly, in all tumor growth predictions the calculated prediction errors (20) were significantly low, as it can be seen from their numerical values in Table V. The tumor growth prediction curves for all the cases examined are shown in Figure 7.
Flow chart of the proposed procedure for the measurement and the TGI mathematical model update for the effective adaptive short term ahead tumor growth prediction.
Average tumor masses registered in the experiment of Ref. (14).
TGI mathematical model best fitting errors.
Adaptive prediction of the tumor growth inhibition. 1: TGI model's parameters estimation until te=25. 2: 1, 2, 3, 4, 5 time periods ahead tumor growth predictions. 3: Moving window's measurements updating and TGI model's parameters re-estimation (i.e. TGI model's identification updating).
The pharmacokinetic data of gemcitabine given as repeated i.p. doses at 100 mg/kg dose level, i.e. the plasma concentration (15).
Discussion
The Tumor Growth Inhibition (TGI) dynamic nonlinear input-output mathematical model introduced by Simeoni et al. in 2004 (10, 11) is a well-established PK/PD nonlinear mathematical model, describing with accuracy the tumor growth inhibition dynamic process in tumor bearing mice and even in clinical practice following administration of chemotherapy drugs. The study in this article was referring first to the effective modeling and numerical calculation of the average tumor mass growth inhibition in mice treated with gemcitabine. The corresponding TGI dynamic model's principal parameters values were estimated effectively by minimizing an objective function formed by the least squares error sum between simulated and experimental TGI values (see Figure 1 and eqs. [2], [3], [4]). A novel and efficient non-linear global optimization technique based on recent deep learning neural networks (DLNN) was successfully used for this purpose (16).
Observed/interpolated (•) and best fitted TGI model tumor growth curves obtained in mice given i.p. doses of gemcitabine (100 mg/kg on Days 19 and 26) using the median (___) and the mean (---) of the estimated values of each unknown parameter.
Prediction errors for the five cases, i.e. one, two, three, four, five time periods (days) ahead, adaptive tumor growth inhibition prediction.
The results showed an almost excellent fit to the experimental data (Figure 6 and Tables III and IV) confirming the ability of the proposed parameters' estimation method to efficiently estimate their values. It was also demonstrated that the median of the repeated estimates performed for each unknown parameter is the best estimate of it.
An innovative adaptive prediction procedure was subsequently described and applied successfully. It is demonstrated, based on the available tumor growth experimental and interpolated data, that it was able to predict with great accuracy the future course of tumor progression for a specific anticancer treatment. This procedure could thus be a significantly useful prognostic tool that could aid towards more personalized and ultimately more successful chemotherapeutic strategies in the clinical practice. Using the presented moving window-like technique on each time instant available experimental data, short term tumor growth prediction for one, two, three, four, five days ahead can be achieved with small and statistically acceptable prediction error. This important outcome is achieved because the TGI mathematical prediction model's principal PD parameters are re-estimated (i.e. updated) at each time instant that new laboratory or clinical practice measurements (observations) are available (i.e. TGI model updating). As was foreseeable, the longer the prediction period of the tumor growth becomes, the lower becomes the convergence of the predicted to the experimental tumor growth curves, reducing the prediction accuracy and increasing the prediction errors. However, in all cases examined, i.e. one, two, three, four, five time periods ahead tumor growth predictions, the square errors observed were relatively small (less than 7.5%) revealing the particularly good short-term prediction capability of the proposed adaptive procedure.
Observed and predicted tumor growth curves obtained in mice given i.p. doses of Gemcitabine (100 mg/kg on Days 19 and 26). One, two, three, four, five time periods ahead adaptive predictions after the 25th day.
Concluding, a novel computational approach is described herein that could be a very helpful tool for predicting the usefulness of a given drug as a putative anticancer agent and to evaluate the expected effect of active doses in humans from preclinical studies using the tumor growth data from human-to-mouse patient derived xenografts (PDX). As a result, several unnecessary studies in mice, costly and time consuming for drug development, could be avoided and the most discriminative dosing regimens, with much lower economic cost, could then be prioritized in much shorter time period. Thus, valuable time for the cancer patients could subsequently be saved and more efficient dose schedules and treatments could be derived, applied, evaluated and updated by performing several preclinical no cost simulations using the presented adaptive TGI predictions. Furthermore, metronomic chemotherapy dosing, which is gaining popularity (21-23), could be investigated and real-time evaluation and adjustment of the antineoplastic agents and their dosing schedules could be further carried out by applying the proposed TGI prediction procedure at each time instant. Obviously, more studies are needed towards this direction but we could anticipate that, ultimately, the choice of the antitumor drugs for each patient could be significantly and more efficiently improved by properly applying the proposed mathematical/computational TGI's prediction approach described herein in the clinical practice.
Acknowledgements
This research has been co-financed by the European Regional Development Fund of the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship and Innovation, under the call RESEARCH–CREATE–INNOVATE (project code: T1EDK-01612).
Footnotes
Authors' Contributions
GSS and KSD conceived the study. SGL, GSS and KSD participated in the design of the study. SGL performed the mathematical analyses. GSS and SGL analyzed and/or interpreted the data. KSD provided the xenograft data. SGL, GSS and KSD wrote the manuscript.
Conflicts of Interest
The Authors declare that they have no conflicts of interest.
- Received July 10, 2020.
- Revision received July 24, 2020.
- Accepted July 27, 2020.
- Copyright© 2020, International Institute of Anticancer Research (Dr. George J. Delinasios), All rights reserved