OPTIMAL CONTROL FOR A BONE METASTASIS WITH RADIOTHERAPY MODEL USING A LINEAR OBJECTIVE FUNCTIONAL

. Radiation is known to cause genetic damage to highly proliferative cells such as cancer cells. However, the radiotherapy eﬀects to bone cells is not completely known. In this work we present a mathematical modeling framework to test hypotheses related to the radiation-induced eﬀects on bone metastasis. Thus, we pose an optimal control problem based on a Komarova model describing the interactions between cancer cells and bone cells at a single site of bone remodeling. The radiotherapy treatment is included in the form of a functional which minimizes the use of radiation using a penalty function. Moreover, we are interested to model the ‘on’ and the ‘oﬀ’ time states of the radiation schedules; so we propose an optimal control problem with a L 1 -type objective functional. Bang-bang or singular arc solutions are the obtained optimal control solutions. We characterize both solutions types and explicitly give necessary optimality conditions for them. We present numerical simulations to analyze the diﬀerent possible radiation eﬀects on the bone and cancer cells. We also evaluate the more signiﬁcant parameters to shift from a bang-bang solution to a singular arc solution and vice versa. Additionally, we study a fractionated radiotherapy model that yields an output solution that resembles intermittent radiotherapy scheduling.


Introduction
Primary tumors may develop secondary metastasis tumors at distant organs, which are the major culprits causing cancer patient death.Frequently, organs like brain, lung, liver and bone are targets of metastatic colonization and in particular, breast and prostate cancers are more prone to develop bone metastasis.These diseases compromise patients' life quality due to severe complications such as fractures and spine cord compression [28,37].These bone-related complications arise because of an intricate interaction network between invasive cancer cells and bone cells.
Metastatic cancer cells interact with the bone microenvironment cells by corrupting a single bone remodeling process [37].More specifically, cancer cells disrupt the communication between osteoclasts and osteoblasts cells, which are in charge of the bone remodeling unit (BRU).Due to this distortion in the communication of the BRU cells, serious bone diseases are developed.For example, bone metastasis breast cancer patients usually present osteolytic lesions, whereas bone metastasis prostate cancer patients show osteoblastic and mixed lesions [14].Therapeutic strategies for bone metastasis include bone-targeted therapies used to diminish excessive osteoclast or osteoblast bone-changing activities.Besides these therapies, cancer-targeted strategies such as chemotherapy and radiotherapy may also be incorporated to control metastatic bone disease [28].
It is known that radiation therapy is used to cause fatal DNA damage to tumor cells.However, DNA damage may occur also to healthy cells which are more able to recover than cancer cells [23,35].Thus, radiation affects osteoblast and osteoclast dynamics [56], but the radiation effects on bone cells are still under constant study since the radiobiological mechanisms are elusive [12,13,52].The most common type of radiation therapy delivered to bone metastasis patients is the external beam radiation therapy [15,17].Several studies have focused on comparing single fraction against multiple fractions in reducing cancer-associated bone pain, and the recommendation favors single fraction for uncomplicated bone metastasis.Despite evidence-based recommendations for radiotherapy scheduling many oncologists favor multiple fraction deliver of radiation [5].Thus, determination of the best radiotherapy scheduling for bone metastasis is still an important issue.
From a general perspective, personalization of radiotherapy scheduling via computational and mathematical optimization methods have been used over the years.There are two well-established basic models to estimate the radiotherapeutic doses.The first one is the linear-quadratic (LQ) model [8], and the second one is referred as the biologically effective dose (BED) calculator [21].The LQ model expresses the surviving cell proportion as a quadratic exponential formula in terms of radiation dose.BED is a useful expression that connects the effectively delivered radiation effect in terms of total dose and a ratio of radiosensitivities parameters.However, some recently works based on differential models studying temporal effects of radiotherapy have been focused their efforts on analyzing several radiotherapy schedules and dose distributions [18,40,44] and searching for optimal treatment schedules [9,11,55].
Another way to deal with the cancer treatment optimization problem is using control theory via differential equations [45,47].Some relevant works based on this focal point are: In Ergun [19] and Ledzewicz et al. [32], they studied optimization of radiotherapy in combination with anti-angiogenic therapy.Rocha et al. [43] proposed an optimization scheduling for a combination of radiotherapy, chemotherapy and immunotherapy.In Ledzewicz et al. [30] the authors analyzed a sub-optimization approach to obtain more clinically relevant chemotherapy treatment delivering.In Ledzewicz and Schättler [29] and de Pillis et al. [16], continuous and switch-like optimal chemotherapies were obtained.In the context of bone metastatic diseases there are few works related to treatment optimization, for example: In Warman et al. [54] based on an evolutionary game model describing prostate cancer bone metastasis progression, they studied optimal chemotherapy scheduling.Lemos et al. [33], Camacho and Jerez [10] and Jerez et al. [26] considered an optimal control problem of the Komarova et al., see [27]; the goal of the first one was to obtain a therapy to accelerate bone mass recovery in the case of the myeloma bone disease, and for the second and third work the goal was to study the efficiency of the radiotherapy, denosumab and antigen receptor T-cell treatments by proposing continuous control solutions in metastatic bone disease.
In this work we present a mathematical model using the optimal control approach to study radiotherapy scheduling optimization in metastatic bone disease.Unlike previous works, we consider a L 1 -type objective functional with the interest to describe the application of intermittent doses of treatment such as in radiation therapy.It is known that the possible solutions for this optimal control problem are the bang-bang and singular arc solutions.A bang-bang solution is a piecewise function which well describes a single dose of radiotherapy.On the other hand, singular arc solutions may not represent realistic doses functions.Thus, we will study the characterization and optimality for these two solution types to find out which are the most relevant model parameters to go from one optimal solution to another.By considering an osteolytic bone lesion, we will analyze the efficiency of the proposed optimal control problem based on minimizing the damage to bone cells and controlling/reducing tumor growth.We will also evaluate numerically some model parameters that are particularly significant to shift from a bang-bang solution to a singular arc one.Finally, we will also propose a simple model based on radiotherapy fractionated scheme to simulate a complete schedule of radiotherapy doses.
The organization of this work is as follows: In Section 2, we present an ordinary differential system with power-law growth functions that models the 'vicious cycle' of bone metastasis with the incorporation of radiotherapy effects.In Section 3, we propose an optimal control problem with a linear term in the objective which gives rise to bang-bang and singular solutions.A theoretical study of the optimality of the solution is carried out.In Section 4, some numerical simulations are performed for the optimal control problem in presence of osteolytic bone metastases where a low-dose drug radiation scheduling is considered.A fractionated radiotherapy simulation is also illustrated.Finally, in Section 5 we present conclusions of this work.

Mathematical model
There are two main paradigms related to bone metastasis.The first one is the so-called 'Seed & Soil' hypothesis that states that metastatic cancer cells have predilection to invade fertile tissues [39].For instance, breast and prostate cancer patients have high chance to develop bone metastatic tumors [28].The second paradigm is the bone metastasis 'vicious cycle' that describes the main steps of bone metastatic colonization: cancer cells express factors that alter homeostatic bone dynamics, and in turn the unbalance of bone cell activities promotes further cancer progression [37].Our mathematical model relies on the vicious cycle paradigm to describe dynamical interactions between cancer cells, osteoclasts and osteoblasts cells.
Osteoclasts and osteoblasts cells, abbreviated as OCs and OBs respectively, regulate bone mass through coordinated bone elimination and bone deposition activities.According to the vicious cycle paradigm, bone metastatic cells (CCs) break this coordination with a potential secondary tumor growth and bone mass getting compromised.Osteolytic lesions are thought to be the result of CCs-promoted over-activation of OCs, thus leading to a rapid bone mass elimination [37].Osteoblastic lesions are characterized by CCs-promoted overactivation of OBs but this cross-talk is largely unknown with the CCs-OCs one [38,51].The mathematical model that we present next codifies these interactions by simplifying the network of chemical reactions using power-law functions (Komarova type models) [10,25].Let C(t), B(t) and M (t) denote the density of OCs, OBs and CCs at time t, respectively.The model proposed is where α i and β i are activities of cell production and removal, respectively; the exponents g i are the net effectiveness of bone cell paracrine factors [27].The parameters σ j are proportional interaction rates between bone cells and cancer cells.Finally, the inherit cancer parameter are α 3 which is the growth rate of CCs and K which is the carrying capacity of CCs.The bone metastasis model (2.1)-(2.3)possesses two biologically relevant equilibrium points: a cancer-free and a cancer-invasion equilibria.
Local stability conditions for such equilibrium points are discussed in Camacho and Jerez [10], Jerez and Camacho [25] to know how the cancer-invasion equilibrium changes locally by variations of certain parameters it is very useful for simulations.

Radiotherapy model
Although radiotherapy has long focused to understand important radiobiological mechanisms that trigger cancer cell death, only recently the radiation effects on other elements of the tumor microenvironment have been studied [3].Thus, it is of interest to consider that radiation may have relevant effects on bone cells as for the cancer cell population [52].Based on works like Swierniak [48], Ghaffari et al. [22] and Rocha et al. [43], we incorporate the radiation effect into the death rates of the cancer cell populations and bone cells.Taking into account such hypotheses, we proposed the following model for radiotherapy [10]: ) ) where u(t) represents the radiation effectiveness required to satisfy 0 ≤ u(t) ≤ u max .The coefficients c 1 and c 2 are rescaling parameters of the radiation effect on the osteoclast and osteoblast populations, respectively.
As mentioned in the introduction, several open problems remains related to radiotherapy optimization, in particular in the context of bone metastasis.In the following, we propose an optimal control problem in order to analyze the radiation therapy using a linear objective functional.Such proposal is based on maximizing the elimination of cancer cells and minimizing the damage to the bone environment.

Optimal control problem
As a first consideration, we aim to minimize the use of radiation on a given time window [0, T ].The undesirable state variable in this case will be the CCs population M .Thus, we penalize positive values of this variable.In Camacho and Jerez [10] we explored an optimal control with a L 2 -type objective functional, obtaining continuous optimal solutions.Here, we propose the following optimal control problem: where J is a L 1 -type functional given by The consideration of L 1 -type functional instead of L 2 -type functional is that we are interested in obtaining bang-bang solutions that resemble 'on'-'off' behavior of treatments.
The non-negative weight parameter w measures the cost balance of using the control u and having the presence of variable M .Here, cost refers not only to the economical burden but also to side-effects.For instance, it is well established that radiation therapy may produce nausea, fatigue and damage some healthy tissues and thus doses should be controlled [4].Moreover, we consider the toxicity of the radiation on the bone cellular microenvironment by introducing the elimination rates c 1 and c 2 for the OCs and OBs populations.
The previous OCP is solved under a predefined admissible control set D and a fixed time interval [a, b].In this work we regard the set D as the set of Lebesgue measurable functions defined on [a, b] [34].In this work we also shall refer to optimal control solutions as OCP solutions.

Pontryagin maximum principle conditions
We employ the Pontryagin Maximum Principle (PMP) to obtain necessary conditions for optimality [7,41].Let z(t) = (C(t), B(t), M (t)) T be the vector of state variables.Observe that model (2.4)-(2.7)can be expressed as where the drift vector field is given by while the control vector field is represented by The Hamiltonian for (OCP) is given by where •, • denotes the usual dot product.According to the PMP, an extremal for (OCP) is a pair (u * , z * ) for which there exists a multiplier λ 0 ≥ 0 and an adjoint vector λ(t) satisfying the following adjoint system: From the PMP, we also have additional conditions: (i) λ 0 and λ(t) do not vanish simultaneously for an optimal control u * and its corresponding optimal state z * for (OCP), see for example [42].(ii) H is constant with minimal value along the optimal solution.
The bang-bang characterization of an optimal solution u * for (OCP) is with the switching function defined by Proposition 3.1.A bang-bang optimal solution, u * (t), of the form (3.6) satisfies that u * (t) = 0 on a small neighborhood of T and H(t) = M (T ).
If ψ ≡ 0 on an open interval I ⊂ [0, T ], then the derivatives of ψ vanish as well, and singular arcs may arise in the optimal control solution.The optimality of singular arcs is subject to the well-known Legendre-Clebsch condition.In the following section, we study such condition for the optimal control problem (3.2)-(3.4).

Legendre-Clebsch condition
The first order Legendre-Clebsch condition for a minimization problem is given by [7,45] Assuming that the control is singular on I, we have ψ ≡ 0 on I. Thus, using (3.8) we have a characterization of the optimal control [45]: (3.9) The functions P and Q will be obtained for the optimal control problem (3.2)-(3.4).
The Lie bracket of two vector fields f 1 , f 2 : R n → R m is defined by where D is the Jacobian and • denotes the usual matrix-vector multiplication.Time derivatives of the switching function, ψ, defined in (3.7) can be expressed using Lie brackets.For that, we first calculate the Jacobians of vector fields f (z) and g(z) given in (3.3) and (3.4), respectively: Then Therefore, the Lie bracket [f, g](z) is given by (3.10) The following propositions can be proved by direct calculations.
where L(z(t)) = wu(t) + M (t) is the Lagrangian of objective functional of (OCP).
Using the previous propositions for our switching function (3.7) along with the Lie bracket (3.10), we compute ψ : then the adjoint variable is the unique solution of the system Proof.The second equation comes from ψ ≡ 0 and (3.7).From ψ ≡ 0 and (3.11), we have the third equation.Finally, we can obtain the first equation from the second one and by considering that the Hamiltonian given in (3.5) satisfies H = M (T ).The second order derivative ψ is obtained as follows [45]: and thus, the functions P (z, λ) and Q(z, λ) are defined for (3.9).
The existence and uniqueness of solutions of any dynamic problem are useful theoretical results in order to guarantee that convergent numerical methods will produce admissible approximations.In our particular case, the existence of optimal control solutions has been proved based on standard theoretical results.For further details the reader is referred to Camacho and Jerez [10], Fleming and Rishel [20].

Numerical simulations
Here we present numerical simulations for the radiotherapy optimal control problem (OCP) in order to study the performance of the normal extremals derived from previous sections.For this, model parameters are adjusted to describe an osteolytic bone metastatic lesion, see Table 1.These parameters were proposed to be consistent with real quantitative values corresponding to OCs and OBs levels [27].An exhaustive calibration analysis on model parameters for an optimal control system as (OCP) was carried out in Camacho and Jerez [10].In a general way to model an osteolytic bone metastatic lesion, we should take σ 1 > 0, σ 3 > 0 and σ 2 < 0; the parameter σ 4 is set to zero as to model a null direct influence from OBs to CCs proliferation.

Optimal control simulations
One practical numerical algorithm to find approximations of optimal control solutions is the so-called forwardbackward sweep method (FBSM) [34].However, it could be difficult to obtain convergence with FBSM when there may be singular arcs.Another approach is to use a direct method.We used BOCOP [6] (http://bocop.org), which uses a solver for large-scale nonlinear programming named IPOPT [53] For the experimental analysis, we are interested to explore different outlines assuming u max = 0.1: -First, we analyze different possible radiation effects to bone cells: no-effect case with c 1 = c 2 = 0, an only-osteoclast effect with c 1 = 0 and c 2 = 0, an only-osteoblast effect with c 1 = 0 and c 2 = 0, and an all-cells effect with c 1 = 0 and c 2 = 0. -Second, we also study the effect of increasing the weight parameter w of the cost functional J.At the top of each column in the next figures, we depict the corresponding value w considered for that case.Namely, we considered w = 1000, w = 2000, and w = 5000.
In the no-effect radiation case c 1 = c 2 = 0 (Fig. 1), we observe that the OCP solutions are in an 'on' state from the beginning until certain time.This translates in delivering radiation as soon as possible.Notice that when tumor proliferation diminishes due to the treatment, the OCs and OBs dynamics are more stable quantitatively.Increasing the value of the weight parameter w decreases the period of radiation application.
In the case of the only-osteoclasts radiation effect, c 1 = 5, c 2 = 0, the OCP solutions differ notoriously from the no-effects case, see (Fig. 2).In this case, the three OCP solutions focus on higher peaks of OCs levels.Due to the coupling between OCs and OBs, we also observe that OBs levels are heavily affected, the osteoblasts wave essentially disappears when radiation decreases OCs levels.However, this is the scenario in which the § cancer shrinks more strongly.It is consistent with the combined treatments currently applied, where radiation is used in conjunction with other treatment that inhibits osteoclasts activation.When the value of w increases the period of radiation application decreases.
Let us turn our attention to Figure 3, where the case of the only-osteoblast radiation effect, c 1 = 0, c 2 = 3, is presented.Compared to the previous figure, we can observe in Figure 3 that singular arcs are present for the OCP solutions.That is, in these cases, the control function u may attain intermediate values between 0 and u max = 0.1.The OCP numerical solutions have an 'off'-singular-'off' structure, meaning that from u(t) = 0 the solution transitions to a singular arc ending with u(t) = 0 again.Notice that these solutions do not realistically describe a radiation therapy.Increasing the value of parameter w decreases both the period and amplitude of radiation application.
We observe that there are cases when singular arcs emerge, and we want to study thoroughly the parameters that may dictate the emergence of this type of solutions.Based on the theoretical and numerical studies carried out so far, we focus on the parameters w and c 2 .Considering c 1 = 0 and u max = 0.1, in Figure 4 we show only the corresponding OCP numerical approximations for w ∈ {100, 1000, 5000} and c 2 ∈ {1, 2.5, 4}.In this figure, we notice that c 2 may be an important parameter related to the appearance of singular arcs, together with a high value of the parameter w.For instance, when c 2 = 1, the numerical solutions only display a bang-bang behavior, whereas for c 2 = 4 there are singular arcs for the three possible choices for the value of w.From numerical simulations we may conjecture that singular arcs solutions arise when the parameter c 2 > 1.That is, the effect of the radiation on osteoblast cells might change the qualitative behavior of the optimal solution.

Fractionated radiotherapy simulations
Solutions with singular arcs may not be easily adapted in practice.For instance, in Alvarez-Arenas et al. [1] sub-optimal solutions for a chemotherapy model are obtained by applying a value threshold to an optimal solution with singular arcs, resulting in a solution resembling a bang-bang solution.In some of our numerical simulations for the OCP solutions we get singular arcs.Thus, based on the radiotherapy model (2.4)-(2.6),we propose a non-optimal approach (NOR model) to describe a complete radiotherapy treatment which is named fractionated radiotherapy scheme.This approach is summarized as follows: Let τ be the time separation between radiation applications and let be the treatment duration per session, then the algorithm consists of alternating on phases (radiation) of length , with off phases (no radiation) of length τ .
Using the fractionated radiotherapy scheme, we present simulations related to three radio schedules for uncomplicated bone metastases [50]: -1 fraction of 8 Gy.-5 fractions of 4 Gy (a total dose of 20 Gy).-10 fractions of 3 Gy (a total dose of 30 Gy).In regards to the treatment duration per session, , some works based on mathematical modeling of cancer radiotherapy use instantaneous or a short radiation duration of the order of minutes [2,19,24].Other authors assume duration of radiotherapy sessions of the order of hours or days [22,43].In this work, we assumed a value of = 30 minutes for simulations.Additionally, we explore different possible radiation effects to bone cells: no-effect case with c 1 = c 2 = 0, an only-osteoclast effect with c 1 = 0 and c 2 = 0, and an only-osteoblast effect with c 1 = 0 and c 2 = 0. Results are shown in Figures 5-7.
When c 1 = c 2 = 0, bone cell dynamics are not directly affected by radiation (Fig. 5).Cancer cell dynamics are slowed due to radiotherapy, but this effect is not significant due to its interaction with bone cells, particularly with osteoclasts.
On the other hand, when c 1 = 1 and c 2 = 0 we assume that osteoclasts are directly affected by radiation.Figure 6 shows that indirectly both osteoblast and cancer cell populations decreases compared with the no control scenario.
Finally, considering c 1 = 0 and c 2 = 1 implies that radiation affects directly osteoblast population.In Figure 7, we observe that cancer cell dynamics are slowed due to radiotherapy similar to Figure 5 (c 1 = c 2 = 0).In this case, the regulatory effect by osteoblasts to osteoclasts is decreased due to lower osteoblast levels.Using the fractionated radiotherapy scheme, we observe that the solution presents the same qualitative behavior than in the case of the solution of the OCP.That is, when radiotherapy affects the population of osteoclasts, cancer is reduced more strongly, otherwise the cancer has a slow decreasing.It is important to emphasize that when the time between doses is reduced, it results in an improvement of the treatment efficacy.On the other hand, it is noteworthy that this second approach considers almost instantaneous effects of radiotherapy, similar to an impulsive therapy with a very short-term effect.

Discussions and conclusions
In this work, we presented an approach to optimize radiotherapy for bone metastasis while taking into account the interaction between bone cells and metastatic cancer cells.The employed approach is based on the optimal control theory (OCP), a framework amply studied in a diversity of biological contexts.The OCP approach presented here considers an L 1 cost functional.This leads to radiotherapy schedules as bang-bang and bangsingular-bang numerical solutions.We observed that the OCP solutions were able to control tumor progression as possible while controlling bone cell dynamics.Also, we detected some numerical OCP solutions with singular arcs that modulated the actual required optimal dose below the maximum permitted dose.However, it may not be easy to adapt singular arc solutions in practice [1].Thus, in this work, we also incorporated simple §  radiotherapy fractionation simulations to show how schemes that resemble clinical practice affect dynamical solutions of our model.
There are two main drawbacks in regards the OCP solutions: robustness and time extension.The robustness refers to a good behavior of the OCP solutions when uncertainties in the state variables were presented.The second drawback refers to the fact that it is not straightforward to obtain a new OCP solution when the final time in the integration interval is increased, since simulations must be performed again.One of the possible issues with bang-bang solutions is that they will be in an 'on' state in a defined sub-interval of time rather than distributed in smaller sub-intervals.This fact has disadvantages when scheduling patients treatments, particularly for the cases of chemotherapy and radiotherapy.It is sometimes preferred that the schedules of these types of treatments allow some periods of time where no treatment is applied to allow the patient to recover from prolonged side-effects.
Taking into consideration the last mentioned issue, in a fractionated radiotherapy (FR) solutions, is easy to extend the final time because the algorithm only requires initial conditions.
A virtue of this approach is that resting days are considered which is clinically relevant for the recovery of patients from the treatment side-effects.
However, an issue of the FR scheme is that it does not take into account information about the time evolution of model dynamics.From numerical simulations FR solutions, in some scenarios, may fail to control effectively tumor growth, particularly when it is assumed a direct radiation effect on osteoblasts.This failure in controlling tumor growth may be due to the assumption of instantaneous effects of radiation in the FR scheme.Further studies are necessary to evaluate how to consider the effects of radiation in cancer mathematical models.
There are still several open questions related to the effects of radiotherapy to non-cancerous cells.So this work offers some research avenues to tackle such challenges.For example, we aim to incorporate the more established LQ model into the bone metastasis dynamical model in a similar way like in Rockne et al. [44] and Ledzewicz et al. [30] rather than the simplistic linear elimination term.Also, it is clinically important to add more constraints to the OCP approach such as a fixed limit of how much radiation could be delivered or how much increase/decrease of bone cells could be tolerable.Finally, it has been reported that radiation therapy may enhance other types of treatment like immunotherapy, which in turn may aid towards controlling metastatic bone diseases [46].