Positive Impulsive Control of Tumor Therapy—A Cyber-Medical Approach

Chemotherapy optimization based on mathematical models is a promising direction of personalized medicine. Personalizing, thus optimizing treatments, may have multiple advantages, from fewer side effects to lower costs. However, personalization is a complicated process in practice. We discuss a mathematical model of tumor growth and therapy optimization algorithms that can be used to personalize therapies. The therapy generation is based on the concept of keeping the drug level over a specified value. A mixed-effect model is used for parametric identification, and the doses are calculated using a two-compartment model for drug pharmacokinetics, and a nonlinear pharmacodynamics and tumor dynamics model. We propose personalized therapy generation algorithms for having a maximal effect and minimal effective doses. We handle inter- and intra-patient variability for the minimal effective dose therapy. Results from mouse experiments for the personalized therapy are discussed and the algorithms are compared to a generic protocol based on overall survival. The experimental results show that the introduced algorithms significantly increased the overall survival of the mice, demonstrating that by control engineering methods an efficient modality of cancer therapy may be possible.

Positive Impulsive Control of Tumor Therapy-A Cyber-Medical Approach I. INTRODUCTION C YBER-MEDICAL systems play an important role in modern medicine, and their importance is growing.The application of STEM in medicine offers prosperous results in medical practice.For example, engineering methods can be used for brain fatigue detection [1], [2], [3], prediction of in-hospital death of trauma [4], skeleton maturity assessment [5], [6], or Parkinson's disease diagnosis [7].Systemtheoretic methods are used in several drug dosing problems, like control of anesthesia [8], or control of blood glucose level with artificial pancreas [9].
System-theoretic methods can also be utilized to optimize drug dosing in cancer therapies.The therapies used in conventional chemotherapy usually have a large resting time, i.e., a long time between the injections and large injected doses [10].They use the maximum tolerable dose (MTD) in order to achieve a maximal effect without killing the patient.Another approach is the low-dose metronomic (LDM) therapy, which applies lower doses with larger frequency.In some cases, this approach was proven more effective against cancer cells, which often become resistant to the treatment [11], [12].LDM therapy can also be cheaper with fewer side effects.We aim to provide algorithms for the mathematical model-based generation of LDM therapies.
Scheduling LDM therapy and providing the required doses is a challenging task.A promising engineering approach is to use a mathematical model describing the effect of the drug on tumor growth and use this model to generate an optimal therapy.There are numerous models in the literature (see [10], [13], [14], [15] and several therapy generation algorithms [10], [16], [17], [18]).A specific characteristic of this physiological problem is that the input is the injection, which is positive, and the system is impulsive.Such systems are rare in engineering practice, and thus handling them requires unconventional solutions [19], [20], [21].Besides therapy generation, the usage of nanorobots in cancer treatments is also exploited by Shi et al. [22], [23], [24], while robotic capsules are used for site-specific drug delivery in [25].
c 2023 The Authors.This work is licensed under a Creative Commons Attribution 4.0 License.
For more information, see https://creativecommons.org/licenses/by/4.0/Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Pérez-García et al. [10] used their model and heuristic scheduling based on in silico tests.Tse et al. considered multidrug optimization problems in [26].Cacace et al. [18] solved the therapy generation as an optimization problem to minimize a quadratic cost function combined with a state observer.They generate antiangiogenic therapy based on the model created by Hahnfeldt et al. [27].The authors of this manuscript created an optimization algorithm for the antiangiogenic therapy based on the Hahnfeldt model in [16], where the aim was to track a reference path that results from maximal dosing with a predefined maximal deviation from that path.The authors used an approach based on pharmacokinetic (Section II-C) and pharmacodynamic (Section II-B) model to optimize chemotherapy in [17] and [28], based on the results of Kusuoka et al. [29].Here, we develop the results further and provide several model-based approaches for the therapy generation problem.
The fundamental component of optimization is the underlying mathematical model.Our algorithm depends on the pharmacokinetics and pharmacodynamics of the drug; however, several strategies discussed here also depend on the tumor dynamics.The literature on tumor models is rich; see the works of [30], [31], and [32]; in our work, we use a tumor model described by ordinary differential equations.
We use a fourth-order model (Section II-A) created to describe measurements from animal experiments first for antiangiogenic therapy using bevacizumab [33], [34], and later for chemotherapy using pegylated liposomal doxorubicin (PLD) in [35] and [36] based on the experiments from [37].The latest model has four state variables, two state variables to describe the living and dead tumor volume dynamics, and two state variables to describe the pharmacokinetics of the drug as a two-compartment model.We formulate the optimization problem as maintaining a predefined drug level during the therapy with the lowest amount of injections, i.e., we carry out optimal impulsive control of the two-compartment pharmacokinetic model and propose three strategies to calculate this predefined drug level.
Impulsive control of compartment systems is not new in the literature; e.g., the work of Pierce and Schumitzky from 1976 [38].Kusuoka et al. [29] created an algorithm in 1981 for finding minimal doses that keep the drug level over a specified lower limit, discussed in Section II-D.Our work is based on this algorithm and contributes by defining the lower limit for the drug level, which we will call minimal inhibitory concentration (MIC) after [39] (see Fig. 1).We test the algorithms in vivo using mouse experiments discussed in Section IV-A.
We propose three strategies to define the lower limit of the drug level: two personalized therapies and one robust therapy created for a specific population.The first strategy is a personalized therapy aiming to maximize the effect of the drug using low dosages and only uses the pharmacokinetic and pharmacodynamic model (Section III-B).The second strategy is a personalized therapy, which finds the minimal effective dosages that ensure that the tumor volume does not increase during the therapy for constant and time-varying (but known) parameters (Section III-C).This strategy requires knowledge of the tumor dynamics as well.The third strategy is the robust version of the second one, which uses the minimal and maximal values of the parameters for a specific population and finds the minimal dosages to ensure that the tumor volume does not increase for the worst-case parameter combination (Section III-D).This strategy is based on interval arithmetics, briefly summarized in Section II-F.One can switch between the strategies based on the tumor response automatically, or the switching can be done manually by an expert (Fig. 1).
We carry out in vivo experiments with 29 mice with breast cancer (Section IV) and show with a log-rank test, that our algorithm significantly increases overall survival compared to a generic therapy.

II. METHODS
The mathematical methods used for the therapy optimization are based on the tumor model discussed in Section II-A.The pharmacodynamics of the model is detailed in Section II-B, which is needed to calculate the drug level that has to be maintained during the therapy.The required dosages are calculated based on the pharmacokinetics of the model discussed in Section II-C.Based on the pharmacokinetic model, the minimal dosages required to keep the drug level over a specified limit are calculated as a solution to an optimization problem discussed in Section II-D.In order to personalize the therapy, identification is carried out with a mixed effect model given in Section II-E.Interval arithmetics is reviewed in Section II-F, which will be used to design robust therapy when only the model parameter ranges are known.

A. Tumor Model
The therapy generation algorithm will be based on a fourthorder model describing tumor dynamics, pharmacodynamics, and pharmacokinetics [35], [36], defined by the equations where  compartment (the blood), and drug level in the peripheral compartment (the tissues), respectively.The unit of the volumes is mm 3 , while the drug levels are in mg•kg −1 , compatible with the unit of the injected doses.The parameters of the model are positive and are listed in Table I.The injections are considered as impulsive effects on the central compartment x 3 , i.e., let the K number of injections take place at time instants t k ≥ 0, k = 0, 1, 2, . . ., K − 1 with t 0 < t 1 < • • • < t K−1 and injected doses u k ≥ 0, then there is a discontinuity of the first kind in x 3 at time t k , such that The output of the system is the total tumor volume since x 1 and x 2 cannot be measured separately in the experiments, only the total tumor volume, which is their sum.

B. Pharmacodynamic Model
The pharmacodynamics of the drug is defined by the Hill function x 3 (ED 50 + x 3 ) −1 in (1) and (2), which is a common function used to model the effect of the drug [8], [40].This function expresses that the effect of the drug is saturated, and after a given limit, increasing the drug level yields a very low increase in the drug effect.The median effective dose parameter ED 50 is the drug concentration where the effect is 50%, i.e., the value of the Hill function is 0.5.
The desired value of the drug level in the central compartment that should be maintained by the therapy discussed in Section II-D depends on the value of ED 50 .An important goal of the therapy optimization will be to keep the drug level in the central compartment over the MIC.
One possible strategy, described in Section III-B, is to define the limit as a constant multiple of the ED 50 parameter, i.e., MIC = κED 50 .If κ is sufficiently large, the value of the Hill function is close to 1, i.e., close to the maximal effect of the drug.For example, for the value of κ = 100, the value of the Hill function is 0.99; thus, the drug has at least 99% effect during the therapy.
Another strategy is to calculate the MIC such that the drug prevents the tumor from growing with the lowest dosages during the therapy, which we discuss in Section III-C for personalized treatment, in the case of intrapatient variability for known model parameters, and as a worst-case treatment for a population to treat interpatient variability in Section III-D.

C. Pharmacokinetic Model
The pharmacokinetics of the model described by ( 3) and ( 4) is a linear time-invariant system.The input of the system is positive and impulsive and has an effect on x 3 as described by (5).Since the drug level in the central compartment (x 3 ) has a direct effect on the tumor as described by the term on the right-hand side of (1), the output y p of the pharmacokinetic model is x 3 .Thus, the state-space representation of the pharmacokinetic model is with u being the sum of impulsive inputs that is written using Dirac-delta distributions δ [20], i.e., at t ≥ 0, where K is the total number of injections, t k ≥ 0, k = 0, 1, . . ., K −1 is the time of injections with doses u k ≥ 0, k = 0, 1, . . ., K − 1, and δ is the Dirac delta distribution, i.e., The output of the pharmacokinetic subsystem ( 7) produced from the impulsive inputs (9) at time t can be written as the sum of impulse responses w of the system as In order to avoid ambiguity, we note that from now on, we use the letter w to denote the impulse response of the pharmacokinetic subsystem (and not the washout rate from Table I, which we will not need in the remaining of this article).The impulse response of the pharmacokinetic subsystem at time t ≥ 0 is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where λ 1 and λ 2 are the eigenvalues of the system matrix A in (7) and can be expressed with the parameters as The pharmacokinetic subsystem is kinetic, thus it is positive [41], [42], [43], [44], which implies that the impulse response of the system is also positive for all t ≥ 0. Theorem 1: The pharmacokinetic subsystem ( 7) is asymptotically stable and nonoscillatory [28].

D. Optimal Impulsive Therapy
One way to optimize impulsive therapy is to find the minimal injection doses in order to maintain a predefined drug level in the patient, i.e., look for the optimal injection doses u k given at time instants t k , k = 0, 1, . . ., K − 1 such that x 3 (t) ≥ MIC for all t ≥ 0. This problem has been addressed by Kusuoka et al. [29], who formulated and solved this optimization problem for compartmental systems with a constraint that the input should also be positive.
Let u = (u 0 , u 2 , . . ., u K−1 ) and 1 = (1, 1, . . ., 1) be a column vector with elements of one and length K, and let be the matrix of impulse responses constructed as where i, j = 1, 2, . . ., K, and t K is a time instant after the last injection, i.e., t K > t K−1 .The purpose of t K is to maintain the predefined drug level for a certain time after the last injection.
The goal is to have y p (t) ≥ MIC for all t ≥ 0, where MIC is the desired lower limit for the drug level, with the constraint that the injections are positive, and the goal of minimizing the total amount of injections.Since the pharmacokinetic subsystem is positive, asymptotically stable, and nonoscillatory, y p (t) ≥ MIC for all t ≥ 0 is equivalent to y p (t k ) ≥ MIC for all t k , k = 0, 1, . . ., K. Thus, the optimization problem can be written as The solution to this optimization problem [29] is This constrained optimization problem minimizes the sum of the drug dosages u (i.e., the cumulative dose) with the constraints such that the drug level in the central compartment is not less than the MIC at the time of the injections (before the drug is injected), and the drug dosages are non-negative.

E. Mixed-Effect Model
The application of the therapy generation algorithms requires the knowledge of the model parameters, i.e., prior parametric identification has to be carried out.If the modeled population has similar intrinsic and extrinsic parameters, then the mixed effect model is an efficient tool for parametric identification [45].In brief, such models assume that actual parameters for a subject-mouse-are random variates from a distribution characterized by their mean and variance (and typically assumed to be normal), and the focus of the estimation lies in these population parameters.These approaches capture the intraindividual correlations of repeated measured data and can be considered as a middle ground between estimating individual parameters and one global parameter set.
Individual variates are obtained similarly to residuals.In fact, a mixed model can be written (with a single level of grouping) as where y i is the response variable, Z i describes the grouping structure, b i represents the so-called random effects, while β is the vector of fixed effects, with X i being the usualfixed effects-design matrix, and ε i being the usual error term.Typically, b i ∼ N (0, ) and ε i ∼ N (0, σ 2 I) is assumed [45].In the present case, y i means the measured tumor volumes, X i and Z i contain the times of the measurements, while β collects the parameters of the differential equations (with b i being their associated random effects).The likelihood implied by the above model can be calculated, and maximum-likelihood estimation is readily possible (although often other variants are used).
A further complication in such parameter identification problems is that the model is specified through a system of differential equations from which no explicit formula can be derived for the response variable.The stochastic approximation expectation-maximization (SAEM) approach can be used, which is an extension of the classic EM algorithm [46], widely used to carry out maximum-likelihood estimation (locally) when the model depends on unobserved parameters, as described by Delyon et al. [47].This is widely used to estimate models described by differential equations with mixed effects [48].
Calculations in the in vivo experiments given in Section IV of this study are carried out under R statistical environment version 4.1.0[49] using the nlmixr package version 2.0.4 [50].

F. Interval Arithmetics
Interval arithmetics is an efficient tool to analyze expressions in which the variables are perturbed, but their lower and upper limits are known.We will use interval arithmetics in Section III-D to calculate the min-max therapy in the presence of parametric perturbations, i.e., calculate the minimal doses that ensure that the drug has inhibiting/killing effect all the time in the worst-case parameter combinations.We will use the notations and definitions from the work of Alefeld and Mayer [51].
Let a denote the lower limit of the variable a, and a denote the upper limit of the variable a, and let [a] denote the closed interval of the possible values of a, i.e., [a] = a, a . ( Let • be a binary operation on intervals, such that • ∈ {+, −, •, /} defined as the set Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where 0 / ∈ [b], if the operation is division.The results of the operations are the intervals In the case if a > 0 and b > 0, we have that The multiplicative inverse is defined as which can be calculated as The standard interval functions are ψ ∈ F = {sin, cos, tan, arctan, exp, ln, abs, sqrt}, which are defined by their range as Functions composed of standard interval functions and binary operations can be used in interval arithmetics based on the definitions.
If a function depends on an interval variable more than once, then we should use interval variables defined as the convex combination of their lower and upper limits [52], i.e., Let ψ be a function of a, then with the lower limit defined as the solution of and the upper limit as the solution of with α constrained to the interval [0, 1].

A. Tumor Dynamics
Consider the first two equations ( 1) and ( 2) of the tumor model and denote the drug effect by E = bx 3 (ED 50 + x 3 ) −1 .These equations describe the tumor dynamics The equilibrium points of this model are the solutions to The solutions are as follows. 1) Steady-state 1 is the origin, which is ideally the goal of the therapy.Steady-state 2 is a nonzero equilibrium that is achieved if the drug effect and the tumor dynamics are balanced.Note that in steady-state 2, the variables can be zero, but in this case, we get steady-state 1.
The Jacobian of ( 30) and ( 31) is with the eigenvalues Since the parameters are positive, we have that λ 2 < 0. In steady-state 2, we have λ 1 = 0, thus the steady-state is on the line Note that we are interested in the case when a − n > 0, i.e., the tumor grows without therapy.
Theorem 2: If E(t) > a − n for all t ≥ 0, then the origin of the model ( 30), ( 31) is asymptotically stable.
Proof: Let E 1 and E 2 be time functions that satisfy that E 2 > E 1 > 0, i.e., E 2 (t) > E 1 (t) > 0 for all t ≥ 0, and let x (1) 1 be the solution to (30) with input E 1 and x (2)  1 be the solution with input E 2 to (30) with initial condition x 1 (0) > 0. Due to the positivity of x 1 and E, we have that ∀t ≥ 0 (37) thus x (1)  1 (t) > x (2)  1 (t) whenever t > 0. Thus, the tumor dynamics model is monotonous in E. Let E 1 > a−n with E 1 ≡ const, and E 2 be an arbitrary function with E 2 (t) > E 1 for all t ≥ 0. Since E 1 > a − n, we have that a − n − E 1 < 0, thus the solution to (30) with initial condition x 1 (0) tends to zero; thus, the origin is asymptotically stable.Due to the monotonicity, x (2)  1 (t) < x (1)  1 (t) for all t > 0, thus the solution with E 2 also converges to zero, since if where the right-hand side tends to zero whenever Thus, if the net effect E is kept over the limit a−n during the whole therapy, then the total tumor volume tends to zero.For impulsive systems with linear pharmacokinetics and nonlinear pharmacodynamics, the effect of the drug can be kept over a specified limit with the method discussed in Section II-D.The following section will provide strategies to define this specific limit.In the original model E = bx 3 (ED 50 + x 3 ) −1 , thus we will provide limit for x 3 denoted as MIC.
Since the maximum of x 3 (ED 50 + x 3 ) −1 is 1, the direct consequence of Theorem 2 follows.
Theorem 3: Suppose that a − n > 0. There exists x 3 such that the origin of the model ( 30), ( 31) is asymptotically stable if and only if b > a − n.
Proof: Suppose that b > a−n.In the original model, we use Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Thus, in the applied tumor model, the effect E is in the interval [0, b), provided that x 3 is always positive.Since lim we can always find a function x 3 which satisfies Theorem 2. Now, suppose that the origin of the model is asymptotically stable.We show indirectly that it implies b > a − n.Suppose that b < a − n, and the origin is asymptotically stable.Thus, we have the limits lim Since the function x 3 (ED 50 + x 3 ) −1 is strictly monotonously increasing and a − n − b < a − n due to the positivity of the parameters, the x 1 solution to (1) with initial condition x 1 (0) is bounded as where e (a−n−b)t tends to infinity due to (42); thus, the origin is unstable, which leads to a contradiction.Due to the nonlinearity of E, the fourth-order model described by ( 1)-( 4) can have rich dynamics; it was shown in [53], that if the therapy is considered to be continuous, i.e., we add u to the state x 3 , and the applied control method is P-type control based on the total tumor volume x 1 + x 2 , which is then the closed-loop system has a nontrivial equilibrium, and the system can have bifurcations with realistic parameter values.Moreover, the qualitative property of the equilibrium is independent of k.

B. Personalized Therapy With Maximal Effect
The critical parameter of the therapy design is the minimal drug level (denoted by MIC) that has to be maintained.If the MIC is known, the algorithm discussed in Section II-D only requires the knowledge of the pharmacokinetic parameters, and can be used to calculate the minimal doses.The value of MIC can be acquired intuitively or based on the pharmacodynamic parameters.For example, choosing MIC = κED 50 with κ sufficiently large ensures high efficiency for the therapy.The choice of κ = 100 results in 99% efficiency based on the pharmacodynamic model.
The advantage of the maximal effect is that it is robust against parametric perturbations due to oversized doses.The other advantage is that we only need to know the pharmacokinetic (PK) and pharmacodynamic (PD) parameters and do not require knowledge about the other tumor model parameters.Thus, the method can also be applied to other tumor models with similar PK and PD models.The disadvantage of the method is the high level of toxicity that may result from the large doses.In order to avoid toxicity, the application of lower doses that still have inhibiting/killing effect may be desirable.

C. Personalized Therapy With Minimal Drug Level
One can also calculate MIC as the minimal drug level that ensures that the tumor volume does not increase between injections.This can be calculated in a personalized fashion if the model parameters of the patient are known.After the MIC is acquired, the algorithm defined in Section III-C is used to generate the required doses.
Theorem 4: The therapy ensures that the tumor does not grow using minimal doses if the drug level is kept over MIC = κED 50 , where Proof: The drug level ensuring that the tumor does not grow can be calculated based on (1), by considering that with the optimized therapy x 3 (t) ≥ MIC for all t ≥ 0, thus the net growth rate is which implies that if we ensure that the net growth rate is smaller or equal to zero at the time of the next injection (prior to the injection), then the net growth rate is smaller than zero between the injections since the PK model is asymptotically stable without oscillatory behavior by Theorem 1 [29].Thus, the tumor volume does not increase if we have From this equation, we can express and calculate the minimum value of κ by dividing the righthand side with ED 50 .This result can be used to generate the minimal dose therapy that ensures that the tumor does not grow during the therapy.However, this approach only works if the parameters do not change during the therapy.Now consider the case of intrapatient variability of the parameters, i.e., when the tumor model parameters change during the therapy.Suppose that the parameters are constant between injections and can only change at the time of injections.Let a i , n i , b i , and ED i 50 denote the value of the parameters a, n, b, and ED 50 prior to the ith injection with i = 0, 1, . . ., K − 1, respectively.Let w i denote the impulse response function of the PK model with the PK model parameter values prior to the ith injection, and x i 3 (0 − ) be the drug level at the time of the ith injection, without the injected dose.Note that we will use t = 0 as the time of the ith injection.
Let d i denote the dose of the ith injection and T i denote the time between the ith injection and (i+1)th injection.The drug level in the central compartment at time T i is composed of the drug level due to the ith injection with a value of d i w i (T i ) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
and the drug level due to the previous injections, denoted by x3 (T i ), which can be calculated as if i > 0, and x3 (T 0 ) = 0. Thus, the drug level in the central compartment at time T i is Theorem 5: The therapy ensures that the tumor does not grow with minimal doses in the case of changing parameters if the drug level is calculated as Proof: Between the ith and (i + 1)th injections, we have to ensure that the net growth rate is not positive, which can be ensured if we have for all i = 0, 1, . . ., K − 1. Substituting (51) yields from which we can express d i to get (52).

D. Min-Max Therapy for Uncertain Parameters
In the case of unknown parametric perturbations, we provide the minimal drug level dosing described in Section III-C, considering the worst-case scenario of parameter changes.This way, the therapy can handle interpatient variabilities in a robust manner.Let a, b, n, and ED 50 be the lower limits, while a, b, n, and ED 50 be the upper limits of the parameters a, b, n, and ED 50 , respectively.Let c, k 1 , and k 2 be the lower limits, while c, k 1 , and k 2 be the upper limits of the parameters c, k 1 , and k 2 , respectively.Let d i denote the drug dose of the ith injection and T i denote the time between the ith injection and (i + 1)th injection.
The time evolution of the drug level is characterized by the impulse response of the PK system, which has to be analyzed with interval analysis to calculate the worst-case value of the impulse response at time T i .Let this worst-case value be denoted by w(T i ), which is the lowest value considering all the parameter combinations.Since the variables appear more than one time in the expression (13), we must define the interval variables as the convex combinations of the lower and upper limits of the parameters, i.e., we have Then, the first eigenvalue of the system matrix is while the second eigenvalue is and the impulse response in the new variables is In the worst-case PK parameter combination, the drug depletes fast thus we need the w(T i ) lower limit for the impulse response at time T i after the ith injection, which can be calculated by solving the minimization problem minimize This problem can be solved using numerical optimization for systems with nonlinear PK equations as well by replacing the impulse response function with the solution of the nonlinear differential equations depending on the uncertain parameters.
The worst-case value of the drug level in the central compartment T i time after the ith injection is Theorem 6: The tumor volume does not grow during the therapy for the worst-case parameter combination if the applied dose is Proof: In order to guarantee that the tumor does not grow between two injections, we must have where we supposed that the current injection takes place at t = 0, and the next injection occurs at t = T i .Then, the value of the left-hand side of (65) is in the interval The therapy ensures that the tumor does not grow in the worst case if Application of interval arithmetics discussed in Section II-F yields that Substituting (62) into (69) yields Rearranging this equation to d i and combining it with the inequality (68) yields the result.

A. Experimental Setup
The results were tested on a clinically relevant, genetically engineered mouse model of breast cancer.In this model, Brca1, a DNA repair gene, and p53, a regulator of cell cycle and genome stability, were knocked out in breast epithelial cells.The resulting mammary tumors highly resemble the Brca1-linked, triple-negative, hereditary breast cancer in humans: the molecular, immunohistochemical, morphological, and genetic characteristics are almost indistinguishable from their human counterpart [54].
Moreover, these tumors respond to chemotherapy very similarly; initial treatment with doxorubicin, docetaxel, or cisplatin significantly reduces tumor size and induces remission.However, long-term therapy often fails due to the emergence of drug resistance [55].Despite we showed that PLD increases relapse-free and overall survival by 6-and 3-fold, respectively, these tumors cannot be cured using conventional chemotherapy regimens [37].Findings obtained by using this model are frequently translated to human cancer clinics due to their similarity to human breast cancer.
The ideal therapy is personalized therapy, which is tailored to the specific patient.This requires knowledge of the patientspecific model parameters.However, the identification of the parameters requires several measurements, which is a large drawback of the method.The measurements can be used to carry out the parametric identification using a mixed-effect model described in Section II-E.The acquired parameters can be used to design optimal therapy using the algorithm described in Section II-D.
There were 29 mice divided into three groups.1) Control group (C): Eight mice were treated with a standard protocol.2) Group S1: Eleven mice, received 4 mg•kg −1 drug when the tumor first reached 200 mm 3 , then the therapy was generated using the methods in Section III after tumor relapse.3) Group S2: Ten mice, received 6 mg•kg −1 drug when the tumor first reached 200 mm 3 , then the therapy was generated using the methods in Section III after tumor relapse.The control group (C) treated with the standard protocol received 6 mg•kg −1 PLD, when the tumor reached 200 mm 3 volume, and the injections were repeated every 10 days if the Fig. 2. Overall survival of the mice (in days).Time is calculated from the second injection, and the survival probability is calculated as the ratio of the living mice.The p value is the significance resulting from the log-rank test of the control group and the groups S1 and S2.
tumor volume was still over 200 mm 3 .Otherwise, the treatment stopped until the tumor trigger, i.e., when the volume reached 200 mm 3 again.This protocol uses large doses (the maximal tolerable dose of PLD is 8 mg•kg −1 ) with relatively large resting times (i.e., not less than ten days).
The mice in the experiments are genetically identical; thus, similar model parameters are expected.They were implanted tumor pieces derived from the same tumor, and their treatment started when their tumor first reached 200 mm 3 .The tumor width and length were measured using calipers three times a week (on Monday, Wednesday, and Friday), and the tumor volume was approximated using the formula [33] V = π 3 (width • length) (2/3) . (71)

B. Results of the Experiment
The tumors reacted well to the first dose (4 mg•kg −1 in the case of group S1, and 6 mg•kg −1 in the case of group S2), and they shrank to a small volume (tumor remission).When the tumors started to grow back (tumor relapse) and reached 200 mm 3 again, we initiated a protocol based on our results presented here in groups S1 and S2.The mice received two injections per week, one on Tuesday and one on Friday.We set 6 mg•kg −1 as the maximal dose that can be given in one single injection.
At the time of the first relapse, we carried out parametric identification for the mice in groups S1 and S2 based on the previous measurements using the algorithm in Section II-E.We generated maximal effect therapies as given in Section III-B using κ = 100, i.e., the minimal drug level was 100 times the median effective dose, resulting in 99% effect of the drug.
We repeated the parametric identification on day 82 of the experiment.Based on expert knowledge, we changed the treatment strategy based on the following heuristics.
1) If the tumor volume was still over 200 mm 3 , we continued applying the maximal effect therapy.If the resulting doses were larger than 6 mg•kg −1 with the new parameters, then we decreased the parameter κ to 10, i.e., the new minimal drug level was ten times the median effective dose, resulting in 90% efficiency.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.2) If the tumor volume was under 200 mm 3 , we switched to minimal drug level therapy described in Section III-C.For comparison, we calculated the survival time of each specimen starting from the second injection since this is the time when the different treatment protocols started.The overall survival is shown in Fig. 2 for the three groups.The survival curves were estimated with the nonparametric method of Kaplan and Meier [56].The equality of survival curves among different groups was tested with the nonparametric log-rank test [57], which resulted in a p-value of 0.031.This analysis showed that the strategies described in this article significantly increased the overall survival of the mice.
We show measurements for four selected mice in Figs.3-6.The figures show the approximated tumor volumes, measured mouse weights, and injected doses.Two mice are  selected from groups S1 and S2, and two from the control group C. Compared to the control group, the injections are more frequent in groups S1 and S2 and have smaller doses.
In these examples, the maximal tumor volumes were larger after relapse for the mice from groups S1 and S2, possibly due to the lower dosages of the drug.However, after the tumor remission, the tumor volume was kept at a low value with small dosages, while in the case of the control group, there were no injections after remission; thus, the tumor grew back after a short time.The results illustrate that the tumor can be kept at a low size with more frequent but smaller doses.
Based on these results, a suggestion for the application of the results is as follows.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
2) Carry out parametric identification using the algorithm presented in Section II-E.3) If the parameters of a specific patient are known, calculate the maximal effect therapy discussed in Section III-B.If toxicity is an issue, either decrease the effect by decreasing κ, or use the minimal dose therapy discussed in Section III-C.4) If the patient-specific parameters are not known, but we have a priori measurements from a population with similar characteristics, we can use the ranges of the parameters to calculate the worst-case therapy discussed in Section III-D.

V. CONCLUSION
Cyber-medical systems in therapy generation may revolutionize the treatment of numerous diseases, including cancer.We proposed therapy generation algorithms based on a tumor model, supposing that the parameters of the model are known from an identification procedure.
The algorithms are based on an optimization that calculates the minimal injection doses, ensuring the drug level is kept over a specified limit based on personalized pharmacokinetic model parameters.We proposed strategies to determine this limit.
The maximal effect therapy is used to keep the drug level over a constant multiple of the median effective dose parameter to ensure that the drug has a predefined efficiency all the time.The advantage of this method is that it only requires extra knowledge of the pharmacodynamics, and by choosing a large multiplier, the robustness can be increased at the cost of toxicity.
We provide a formula for the minimal drug level that ensures that the tumor volume does not increase during the therapy.These results were generalized for the cases when intrapatient and interpatient variability of the parameters are present, providing a robust impulsive therapy solution.
The results were tested in in vivo experiments using 29 mice.The therapy generated with the proposed methods was compared to a generic protocol used in clinical practice.The superiority of our results was proved with a log-rank test, which showed a significant increase in overall survival.
The combination of the strategies given in Section III was carried out using heuristics in the experiments.This can be developed further by using feedback of the tumor volume, which can change the MIC after every sampling, or one can track the parameter changes and use the strategy in Section III-C to personalize the therapy and track intrapatient variability.

Fig. 1 .
Fig. 1.Procedure of the therapy optimization: the MIC is specified based on a chosen strategy, and the minimal doses are calculated, which ensure that the drug level is over MIC.

Fig. 3 .Fig. 4 .
Fig.3.Tumor volume measurements, weight measurements, and the injected doses for mouse 6 from group S1.The measurements are denoted with x-es, which are linearly interpolated for visualization.

Fig. 5 .
Fig. 5. Tumor volume measurements, weight measurements, and the injected doses for mouse 1 from group C. The measurements are denoted with x-es, which are linearly interpolated for visualization.

Fig. 6 .
Fig.6.Tumor volume measurements, weight measurements, and the injected doses for the mouse 5 from group C. The measurements are denoted with x-es, which are linearly interpolated for visualization.
x 1 , x 2 , x 3 , and x 4 are the time functions of the living tumor volume, dead tumor volume, drug level in the central Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.