Pharmacokinetics and Pharmacodynamics Models of Tumor Growth and Anticancer E ects in Discrete Time

Abstract: We study the h-discrete and h-discrete fractional representation of a pharmacokineticspharmacodynamics (PK-PD) model describing tumor growth and anticancer e ects in continuous time considering a time scale hN0, where h > 0. Since the measurements of the drug concentration in plasma were taken hourly, we consider h = 1/24 and obtain themodel in discrete time (i.e. hourly). We then continue with fractionalizing the h-discrete nabla operator in the h-discrete model to obtain the model as a system of nabla h-fractional di erence equations. In order to solve the fractional h-discrete system analytically we state and prove some theorems in the theory of discrete fractional calculus. After estimating and getting con dence intervals of the model parameters, we compare residual squared sum values of the models in one table. Our study shows that the new introduced models provide tting as good as the existing models in continuous time.


Introduction
In a typical pharmacokinetics and pharmacodynamics (PK-PD) model, which describes the impact of anticancer treatment on the dynamics of tumor growth, a transit compartment model consists of a system of rst order di erential equations. A rst order di erential equation is an equation that includes an unknown function of time t (call f ) and its rst order derivative, which stands for the instantaneous rate of change of f at time t. In this regard, the PK-PD model can be viewed as a relation between unknown functions of time, their rst order derivatives that include some estimated parameters, and drug concentration as a function of time.
In this paper, we introduce the PK-PD model as a system of h-fractional di erence equations of noninteger order α, where α is any real number between zero and one. A fractional di erence equation is an equation which includes an unknown function of time (call f ); and is de ned as the rate of change of the ( − α)-th order sum of the function f at time t. From a practical point of view, the mentioned non-integer order sum can be considered as the history of the function f from the initial time t to t. Hence it is worth mentioning that the fractional di erence of a function depends on function's whole time history from initial time t to t; and not on its instantaneous behavior at time t. Additionally, the models with fractional di erence equations in discrete time provide the possibility of expressing the solution as an iteration schema. Hence there will be no need to search for numerical approximation methods to solve the equations. Because of such characteristics, we believe that introducing a system of fractional di erence equations into the PK-PD research area will serve as a base for more applicable models in cancer research. In contrast to the previouslyintroduced PK-PD models with ordinary di erential equations (ODE) and fractional di erential equations (FDE) in the literature [1,2,8,9,10,12,13,14,15,17,18,19,22,23,24,25,26,27,28,29], here we take a di erent approach which yields new insights in modeling.
In the paper [5], the authors introduced and then fractionalized the following model on the time scale where u(t) = ln x (t). While the models in [5] were serving as good candidates for the purpose of data tting, some limitations on the models were pointed out such as the discrete model (h = ) looses some necessary properties that hold in the continuous model, a jump occurs in the h-discrete model (h > ) after drug administered. In this work, we aim to keep the data tting capability of the models and to remove the mentioned limitations by introducing and studying the following model in discrete time (h > ): where u(t) = ln x (t). 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 h-discrete fractional equations. The nabla h-fractional di erence operator is considered in the sense of Riemann-Liouville de nition of the fractional derivative. In Section 3, we rst demonstrate the construction principle of our continuous tumor growth model and present important model properties. Then we obtain both h-discrete and h-discrete fractional models after discretizing the PK-PD model in continuous time. In Section 4, for estimating parameters, we use Nelder-Mead algorithms in R and MatLab after solving each h-discrete and fractional h-discrete models explicitly. We close this section by presenting our ndings with con dence intervals on the model parameters. In Section 5, we discuss advantages of the models in discrete time.

Preliminaries
De nition 2.1. Let a ∈ R and h ∈ R + . The backward h-di erence operator for a function f : hNa −→ R is de ned by De nition 2.2. [21] For any t, r ∈ R and h > , the h-rising factorial function is de ned by where the quotient is well-de ned. Here Γ(·) denotes the Euler gamma function.
De nition 2.3. [7] Let α > and a be two real numbers. For a function f : hNa −→ R, the nabla h-fractional sum with order α is de ned by We note that if h = , then we have the α-th order fractional sum operator [4] In the rest of the paper, we write ∇ −α a for ∇ −α ,a .

De nition 2.4. [7]
The nabla h-fractional di erence of order α in the sense of Riemann-Liouville is de ned by where a, α ∈ R, n − < α < n, and n is a positive integer.
where u ∈ N a/h . We use De nition 2.4 and i) to obtain the desired identity. Indeed, we have where u ∈ N a/h+ .
The following theorem states the variation of constant formula in hNa. This theorem has been stated and proved for any time scale in [11].

and t is any real number. Then the rst order nabla h-di erence
has the general solution where c is constant.
The following theorem states the variation of constant formula for α-th order discrete fractional equation in Na. For the proof of the theorem we refer the reader to an article by Atıcı and Eloe [3]. Let α be a real number between zero and one.
Theorem 2.6. Let A be a n × n constant matrix with eigenvalues of modulus less than one and suppose f is a vector valued function. Then the initial value problem has a unique solution. Moreover, this solution is given by Now we are in a position to state and prove the variation of constant formula for fractional h-di erence equations.
has the general solution where c is constant and Proof. Replacing t by uh and using Lemma 2.1-(ii), we write the Equation (2.5) as a fractional di erence equation for u ∈ N h . Then we use Theorem 2.6 to obtain the desired result.
Next we consider the following initial value problem where |λh α | < . By Theorem 2.7, the unique solution of this IVP is for t ∈ hN .

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 [22] 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 [23]. 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 with h > . Hence we have where u(t) = ln x (t). We solve the above system of di erence equations by using Theorem 2.5. Hence we have The perturbed h-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 α ∈ ( , ). We use Theorem 2.7 as a tool to obtain the following solutions.

. The construction principle for the discrete equations
We brie y explain the idea behind the construction principle for the discrete equations we introduced here and in the paper [5]. We consider two possible discretization techniques: T1. First order derivatives in the continuous model can be approximated by use of the ∆ h operator, i.e + h), the discrete equations can be written with ∇ h x(t) following by h unit back shift on the time dependent variables on the right side of the equations.
T2. First order derivatives in the continuous model can be approximated by use of the ∇ h operator, i.e In [5], the authors employed the technique T1 to discretize the rst order di erential equations in the continuous model. Later in the simulations, some limitations occurred in the discrete models. One of the concern was the jump in the hourly model when the drug concentration c(t) was multiplied by some numbers 10, 20, 100, as it was illustrated in Figure 3 in [5]. In this paper, we employ the technique T2. In the later sections, we justify that the new discrete models are as good as the continuous model in data tting, and even better for some cases. As we illustrate in Figure 1, the new discrete models address all shortcomings listed in [5]. In this sense, the discrete models which are introduced in this paper can be considered as the improved PK-PD models in discrete time.

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

Parameter estimates
We assume that n scalar longitudinal observations are represented by the statistical model 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 represents 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 from the tumor growth inhibition model, 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 Since the realized data for both unperturbed and perturbed tumor growth are available, both sets are utilized simultaneously in tting the growth curves (with and without drugs) and obtaining estimated parameters. The unperturbed growth curves are similar to the perturbed counterpart excepts for concentration c(t) = for all t.
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 bŷ 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 based on the exact solutions for the h-discrete and h-discrete fractional models. The estimates of parameters can be obtained by minimizing the residual sum of squares SS using the statistical software package R. Then the sensitivity matrix χ(β) can be derived. The point estimates and the 90% con dence intervals of the parameters for the discrete model are included in Table 1 and those for the discrete fractional model are in Table 2. Table 3 gives residual sum of squares RSS values for all drugs being tted using continuous, discrete, and discrete fractional models using unperturbed and perturbed data simultaneously.
Comparing the RSS values between the discrete model and the discrete fractional model for the same drugs, we can see that the discrete fractional models give better t compared to the discrete counterpart for three drugs A2-120, B-100, and C-150 . In three out of ve drugs, the discrete models gives slightly better t than the continuous models. Note that the estimated parameters and thus RSS for the continuous model are obtained from the numerical ODE solver, while the estimated parameters and RSS are obtained from tting the exact solutions for the discrete and discrete fractional models.

Conclusions
Calculus on a time scale hN serves as a bridge between calculus on R + and calculus on N . In applied mathematics, there are many problems where h-di erence equations are obtained by approximating either ordinary di erential equations (ODEs) or partial di erential equations (PDEs). The main goal is to simplify the problem so that the qualitative behavior of solutions of ODEs or PDEs can be obtained accordingly.
In this study, we discretize the PK-PD model in continuous time. Since the concentration of each drug was measured hourly, we have hN as the domain of the discrete model. We introduce a new set of PK-PD models which behave similar to the PK-PD models in continuous time. These new models bring some advantages compared to previously known models. They are elementary to work with and they also are quite accurate in data tting. The latest techniques in data tting require us to use ODE solver which applies some numerical approximations to the model in continuous time. Here the new models do not require us to use any approximation techniques. We solve them analytically and use their explicit solutions directly in our codes. As we demonstrated in Figure 1 the new models posses the fundamental properties of PK-PD model when we increase the dose.