A study on discrete and discrete fractional pharmacokinetics-pharmacodynamics models for tumor growth and anti-cancer e ects

Abstract: We study the discrete and discrete fractional representation of a pharmacokinetics pharmacodynamics (PK-PD)model describing tumor growth and anti-cancer e ects in continuous time considering a time scale hNh 0, where h > 0. Since the measurements of the tumor volume in mice were taken daily, we consider h = 1 and obtain the model in discrete time (i.e. daily). We then continue with fractionalizing the discrete nabla operator to obtain the model as a system of nabla fractional di erence equations. The nabla fractional di erence operator is considered in the sense of Riemann-Liouville de nition of the fractional derivative. In order to solve the fractional discrete system analytically we state and prove some theorems in the theory of discrete fractional calculus. For the data tting purpose, we use a new developed method which is known as an improved version of the partial sum method to estimate the parameters for discrete and discrete fractional models. Sensitivity analysis is conducted to incorporate uncertainty/noise into the model. We employ both frequentist approach and Bayesian method to construct 90 percent con dence intervals for the parameters. Lastly, for the purpose of practicality, we test the discrete models for their e ciency and illustrate their current limitations for application.


Introduction
Mathematical modeling of pharmacokinetics (PK) and pharmacodynamics (PD) deals with (i) the distribution and elimination of a drug, the pharmacokinetics [17], and (ii) the therapeutic e ect of a drug on a speci c target, the pharmacodynamics [24]. Mathematical PK-PD models became an essential part in drug development and clinical pharmacology in the last 20-30 years. The U.S. Food and Drug Administration recognized such computational modeling and simulation tools as an improvement of the e ciency for developing safe and e ective drugs [18].
Historically, clinicians and pharmacists [15], [16] realized very early that mathematical models are useful tools to answer pharmacological questions. However, it is important to point out that the development of the mathematical PK-PD modeling theory was quite independent from the mathematical community. PK-PD modeling is a eld where the representation of pharmacological concepts in the models is the main focus. Therefore, clinical interpretation of model parameters can be made and reasonable model behavior can be used for simulations. In the last 2-3 decades, mathematicians started to contribute with valuable input to the PK-PD eld. Moreover, a small mathematical community that promotes more advanced mathematical and computational approaches has been established within the PK-PD eld [25].
The main purpose of a PK-PD model is to be a tool to answer speci c clinical questions. More precisely, Allen et al. [1] pointed out recently that the success of mathematical modeling, for example in the biopharma industry, depends on getting the (i) right question, (ii) right model, and (iii) right analysis. Usually, PK-PD models are formulated by ordinary di erential equations and can be developed in a more data-driven approach, i.e. simply characterizing available data, or in a more mechanistic approach, where detailed underlying physiological mechanisms are represented. In this paper, we consider a PK-PD model for tumor growth and anti-cancer drug e ects that can be considered as a hybrid of both approaches, which is often also called a semi-mechanistic model. In general, the complexity of the PK-PD model structure strongly depends on available data, prior knowledge and most important on the speci c questions that need to be answered.
Mathematical modeling of tumor growth started in the 1960s by A. K. Laird [23] with the application of the Gompertz function [19] to t data from di erent animals, see e.g. Bonate et al. [13] for more details. It turned out that this sigmoidal curve ts well to the general tumor growth behavior which usually has three phases, exponential growth in the beginning followed by linear growth and nally saturation. Nowadays the Gompertz function is a common tool in PK-PD modeling although it is a more data-driven than mechanistic approach. The major goal in anti-cancer drug development is to assess the anti-cancer e ect of di erent drug candidates and consequently to early identify promising drug candidates Therefore, the PK-PD models have to characterize the drug potency and need the capability to perform simulations of tumor growth behavior for di erent doses. For these reasons, tumor growth models were extended with (i) PK models describing the drug concentration, and (ii) e ect models relating the action of the drug on the tumor growth [26], [21], [22] and again see Bonate at el. [13] for a broader overview.
Development of advanced and e ective mathematical models which describe and predict tumor growth kinetics is an ongoing research area for scientists, especially for applied mathematicians [2], [6]- [10], [14]- [19], and [21]- [29]. Motivated by recently studied models in fractional calculus [14,28], we aim (i) to discretize a PK-PD model in continuous time and obtain a system of fractional di erence equations in discrete time, and (ii) to investigate the properties of the discrete and discrete fractional model representation.
For the purpose of analyzing the estimated parameters, sensitivity analysis is conducted to incorporate uncertainty/noise into the model. Uncertainty could be resulted from measurements and/or intrinsic variability in experimental units and environment. By accounting for uncertainty in a mathematical model, applied mathematicians can attach measures of reliability to estimated quantities using experiment data sets. The standard procedure involves obtaining con dence intervals for model parameters. In this study, we employ both frequentist approach and Bayesian method to construct 90 percent con dence intervals for the parameters.
The plan of the paper is as follows: In Section 2, we give some preliminaries so that the reader will be familiar with the mathematical formulations in the later sections. We state and prove some theorems that will serve as main tools to obtain the explicit solutions of the system of discrete fractional equations. In Section 3, we rst demonstrate the construction principle of our continuous tumor growth model and present important model properties. Then we obtain both discrete and discrete fractional models after discretizing the PK-PD model in continuous time. For estimating parameters, we use the improved version of the partial sum method after solving each model explicitly. A pseudo-code for the discrete fractional model outlines the algorithm we write in Mathematica-11. We close this section by presenting our ndings with sensitivity analysis on the model parameters. In Section 4, we discuss limitations and advantages of the discrete model representation.
In the models, we consider a nabla (backward) discrete fractional operator of Riemann-Liouville type. We refer the reader to some recent work by Atıcı and co-authors [3]- [5] for the use of the discrete fractional operators in modeling and a book by Goodrich and Peterson [20] for the theory of discrete fractional calculus.

Preliminaries
In this section, we present su cient fundamental de nitions and formulas so that the article is selfcontained. In addition, we state and prove some theorems in discrete fractional calculus so that we explicitly solve the discrete fractional models introduced in the later sections.

De nition 2.2.
Let µ be any real number. The rising factorial power t µ (read 't to the µ rising') is de ned as where t ∈ R \ {..., − , − , }, µ = and Γ denotes the Gamma function. We consider the ν-th order fractional sum of f de ned as where ν ≥ , ρ(t) = t − , and t ∈ Na. Further, we consider the ν-th order fractional di erence (a Riemann-Liouville fractional di erence) of f de ned by where ν > , n − < ν < n, n denotes a positive integer. Let λ, c ∈ R and ν ∈ ( , ). We study the following the initial value problem (IVP) Proof. We use the de nition of the fractional nabla di erence operator (2.2) to obtain the following iteration schema.

Theorem 2.4. The unique solution of the initial value problem (2.3)-(2.4) is given by
Next we interchange the order of sums and obtain By using the formula t r = Γ(t + ) Γ(r + )Γ(t − r + ) and the de nition of the rising factorial power we get, where we used the identity Next we apply the following rule to the above expression Hence, we have = cλy(t − ). Uniqueness of this solution follows from Theorem 2.3.
has the general solution where c is constant and We show that Using the de nition of the nabla fractional di erence operator we have Next we interchange the order of sums and obtain Next we apply the following rule to the above expression Hence, we have Thus, We use Theorem 2.4 to complete the proof.
Theorem 2.6. Assume λ ∈ R\{− }. The rst order nabla di erence equation has the general solution where c is constant.
De nition 2.7. Let a ∈ R and h ∈ R + . The backward h-di erence operator for a function f : hN h a −→ R is de ned by Theorem 2.8. [12] Assume λ ∈ R\{− h }. The rst order nabla h-di erence equation ., (2.11) has the general solution where c is constant.
For further reading on time scale calculus which includes hN h a , we refer the reader to a book by Bohner and Peterson [11].

The tumor growth inhibition model with drug e ect . Development of the continuous tumor growth inhibition model
Replacing the unperturbed growth component of the PK-PD model in [21] with the Gompertz growth component, results in the following representation of the PK-PD model in continuous time where a, b, k , k are model parameters to estimate and c(t) represents the drug concentration in plasma described by mono-or bi-exponential PK models [22]. During anticancer treatment it is assumed that the growth dynamics of the tumor will be perturbed by the anticancer drug e ect described with the model parameter k . Due to drug action, proliferating cells become non-proliferating depending on the drug concentration. The model assumes that cells a ected by drug action immediately stop proliferating and pass through apoptotic stages (x , x ) with a rate k before they die. Since these non-proliferating cells still add to total tumor mass, total tumor volume w(t) is the sum of proliferating x and non-proliferating tumor cells (x , x ).
The tumor growth inhibition model was constructed in such a way that two fundamental properties hold: (P.1) During drug administration, i.e. c(t) > , the tumor growth will be inhibited. (P. 2) The tumor volume will never become negative, i.e. w(t) > for all t ≥ .
We rst discretize the above model considering the time scale hN h with h > . Hence we have where u(t) = ln x (t). We solve the above system of di erence equations by using Theorem 2.8. Hence we have k c(sh)), t = , h, h, ..., If h = , the perturbed discrete tumor growth model with three compartments reads: where u(t) = ln x (t) and a, b, k , k are model parameters to estimate. We solve the above system of di erence equations by using Theorem 2.6. Hence we have The perturbed discrete fractional tumor growth model with three compartments reads: where u(t) = ln x (t) and ν, a, b, k , k are model parameters to estimate. Here we assume that ν is a real number such that ν ∈ ( , ). For the rst equation in the above system, we use some properties of the nabla-di erence operator and the Gamma function, we obtain the following iteration formula for u(t).

. Tumor growth and drug concentration data
Data was taken from the supplemental material of [21]. These data describe the tumor volume measurements in Xenograft mice treated with di erent compounds. For more details see [21].

. Parameter Estimations
In our study we used the improved partial sum method in Mathematica to estimate the parameters for the discrete and discrete fractional models. We include a pseudo-code of the algorithm for discrete fractional model in Algorithm 3 (see Appendix). This algorithm outlines how the improved partial sum method ( [7]) works for the system of equations. In Table 1, we list the parameter estimates and residual sum of squares (RSS) values for the discrete and the discrete fractional PK-PD models. Comparing the RSS values between the discrete and the discrete fractional model, the discrete fractional model results in better tting for all drugs.

. Analyzing the Model Statistically
After the model parameters are estimated, sensitivity analysis is conducted to incorporate uncertainty/noise into the model. Uncertainty could be resulted from measurements and/or intrinsic variability in experimental units and environment. By accounting for uncertainty in mathematical model, applied mathematicians can attach measures of reliability to estimated quantities using experiment data sets. The standard procedure involves obtaining point estimates for the model parameters and then either obtain the corresponding stan-dard error of the estimates or construct con dence intervals. In this study, we employ both the frequentist approach and Bayesian method to construct 90% con dence intervals for the model parameters.
We assume that n scalar longitudinal observations are represented by the statistical model W j = f j (β) + ϵ j ; j = , ..., n, (3.1) where f j (β) is the model for the observations in terms of the state variables and β ∈ R p is a vector of theoretical true parameter values. The error terms ϵ j , j = , ..., n represent noise, measurement error, or uncontrol variables that can potentially in uence the deterministic relationship represented by f j (β). For our statistical model of the observation or measurement process (3.1) we assume that the errors ϵ j , j = , ...n, are independent identically distributed random variables with mean E[ϵ j ] = and constant variance Var[ϵ j ] = σ , where σ is unknown.
Thus we use the realized data w j from observations of three-compartment model to seek a valueβ that minimizes The estimatorβ(W) is a random variable whose realized value depends on the observed data w j of random variables W j . The distribution ofβ(W) (through repeated sampling) is called the sampling distribution. Knowledge of this sampling distribution provides uncertainty information (e.g., standard errors) for the numerical values ofβ obtained using a speci c data set. Under the regularity assumptions, as the sample size n approaches in nity, the sampling distribution ofβ(W) is approximately is a n × p sensitivity matrix with elements χ jk (β) = ∂f j (β) ∂β k ; j = , ..., n; k = , ..., p.
Since β and σ are not known, we must approximate them using the parameter estimates to obtain the estimate for the variance-covariance matrix whereσ is given by Standard errors of the estimateβ k are thus given by SE(β k ) = Σ kk (β), k = , ..., p. Then a ( − γ) % con dence interval for β k is readily given byβ k ±t n−p,γ/ *SE(β k ), where the critical value t n−p,γ/ is computed from the Student's t distribution with n − p degrees of freedom. The 90% con dence intervals for the parameters in the model is computed as follows: The least squares estimatesâ,b,k , andk for parameters a, b, k and k of the discrete model and estimatesâ,b,ν,k , and k for parameters a, b, ν, k and k of the discrete fractional model are obtained from iterative computational scheme. Then the sensitivity matrix χ(β) can be derived. The results of 90% con dence intervals of parameters for the discrete model are included in Table 2 and those for the discrete fractional model are in Table 3.
After obtaining the sensitivity matrix, we utilize the Markov Chain Monte Carlo (MCMC) technique which constructs the Markov chain whose stationary distribution is the posterior distribution of the model parameters. The speci c algorithm is called the Delayed Rejection Adaptive Metropolis [27]. The algorithm starts with an initial vector β and is iterated M times. In each iteration, a new candidate β * is chosen from a Normal distribution that centers at the value β from the previous iteration and the variance-covariance matrix Σ(β) whereβ's are the estimates of the parameters. Whether this candidate is either accepted as the new value for β or not depends on how more likely the candidate is compared to the previous value of β with respect to the posterior distribution. If the candidate is rejected, the delayed rejection algorithm provides a mechanism for constructing alternative candidate β rather than retaining the previous value. Also, during a nonadaptive period of length k , chain values β are computed using the initial covariance matrix Σ(β). Once adaptation starts, the variance-covariance matrix is adapted according to the covariance matrix of the accepted candidates βs in the previous steps. The results of 90% Bayesian interval estimates of parameters for two models are included in Table 2 and 3. Based on Table 2 and 3, we can see that in most cases, the model parameters are signi cant (con dence intervals do not contain zero). Also, the Bayesian approach gives remarkably narrower con dence intervals in most cases considered. The reason is that the construction of con dence intervals using frequentist approach is based upon the asymptotic theory of estimators which assumes that the sample size has to be su ciently large. When the sample size is not large enough, the frequentist asymptotic con dence interval is relatively wide. On the other hand, in the Bayesian inference, the sample size is often not of concern.

Limitations in Discrete Models
Fractionalization of mathematical models in the PKPD eld is not new [28,30] but has not really prevailed. One reason might be the raising mathematical complexity and consequently the di culties in numerical implementation for the user and probably also in run-time. We have to keep in mind that, e.g. in a clinical setting with a population consisting of several hundreds of patients run-time is an issue. As we are presenting here, rst discretizing and then fractionalizing a PK-PD model is an approach which can be resulted in a simple explicit solution for which direct implementation without any numerical solvers is possible. However, simulations of the discrete model (h = ) showed some limitations of the discretiziation approach which will be discussed in this paragraph: First, property (P.1) (tumor growth is always inhibited during treatment) partially fails in discrete models. Simulations showed that for increasing doses the tumor volume can be slightly increased in the very beginning of the treatment and then decreases during treatment as illustrated in Figure 1. In discrete models, because of the nature of the solutions, the e ect of the drug occurs in the following day after the drug is administered. The graph in discrete time does not decrease immediately after drug administration.
Second, property (P.2) (tumor volume will never become negative) holds for the discrete model on hN h with h = , and does not hold for the discrete model (h = ). As in Figure 2, the simulations showed that for strongly increasing doses the tumor volume can become negative. On one hand this is obviously not reasonable from a biological perspective. On the other hand, the model on hN h recovers the property (P.2) as illustrated in Figure 3. In addition to all these, we also want to point out that in practice model predictions are usually only reliable for a local range of parameters.
Third, since the PK is on a di erent time scale than the PD, this may have an impact on the chosen discretization step size. Please note, an absorption compartment in a PK model is necessary, if it takes some time for the drug to reach the blood. Typical examples are orally administered drugs with a delay between administration and appearance in blood of e.g. several hours. In contrast, intravenously administered drugs appear immediately in blood. In this paper, we simpli ed the PK from [21] and omitted the absorption component for two reasons: First, absorption was rapid and can be neglected from a pharmacological perspective. Second, we have to be aware that the drug was administered every day and at the same time when the tumor was measured. Because c(t) describes the drug concentration in blood, c would be zero at the time of administration, since it takes a very short time to reach the blood due to the absorption component. Hence, without an absorption component, c at the time of administration characterizes always the peak of the drug concentration. Because the discrete models always evaluate c every day, this ensures that we always use the maximal drug concentration at the discretization steps. Without this modi cation, c would be always zero in the discrete models. However, it is possible to take such characteristic of c(t) into account in the discrete model on hN h by choosing h in a smaller value.
Forth, the solution representation introduces the restriction of λ ≠ − (as stated in Theorem 2.6) which is then translated to our discrete tumor growth model for the parameters b and k .
In summary, the discrete model (h = ) looses some necessary properties that hold in the continuous model representation. As we illustrated in Figure 3, the second property can be recovered by a model on hN h which reads c(t) in every hour not only one time in a day as in the discrete model. However the graph of the model in hN increases immediately at the time of drug administration starts. This issue may lead to the question that maybe the discretization of the continuous model has to be constructed in a di erent way to assure that (P.1) and (P.2) become valid in the discrete model (h = ) and no jump occurs for the hN h model.   which CT is 0.

Output: a and b
We refer the improved partial sum algorithm