MATHEMATICAL ANALYSIS OF A GENERALISED MODEL OF CHEMOTHERAPY FOR LOW GRADE GLIOMAS

. We study mathematical properties of a model describing growth of primary brain tumours called low-grade gliomas (LGGs) and their response to chemotherapy. The motivation for considering this particular type of cancer is its large impact on society. LGGs aﬀect mainly young adults and eventually result in death, despite the tumour growth rate being slow. The model studied consists of two non-autonomous ordinary diﬀerential equations and is a gen- eralised version of the model proposed by Bogda´nska et al. (Math. Biosci. 2017). We discuss the stability of stationary states, prove global stability of tumour-free steady state and, in some cases, justify the existence of periodic solutions. Assuming that chemotherapy eﬀectiveness remains constant in time, we provide analytical estimates and calculate minimal doses of the drug that should eliminate the tumour for particular patients with LGGs.


1.
Introduction. Low grade gliomas (LGGs), grade II gliomas according to the World Health Organization classification [25], are primary tumours originating from glial cells in brain. They affect mainly young adults and often remain asymptomatic for many years, mostly due to their low growth rate. However, they inevitably transform into more aggressive forms -high grade gliomas (HGGs), causing the neurological deficits and eventually premature patient death for virtually all LGGs patients. Median survival of LGG patients ranges from 5 to 10 years [16], while for HGG patients it is only 1 to 2 years [26].
Typical treatment of LGGs patients involves tumour excision, chemotherapy and/or radiotherapy depending on many factors such as age, performance status, location of tumour and patient preference. It is recommended to perform surgery as extensive as possible due to its evident effect on the patient prognosis [16]. However, the impact and guidelines for the use of chemotherapy are less clear. Chemotherapy tends to be performed early, especially in case of clinical and/or radiological progression. Radiation therapy causes significant neurotoxicity and is proposed usually for inoperable tumours for which chemotherapy proved ineffective and/or at the time of malignant transformation i.e. a transformation of the tumour into more aggressive anaplastic form [42].
There are two effective, approved for clinical use, chemotherapeutic drug/drug combinations to treat glioma patients: temozolomide (TMZ) and a combination of procarbazine, lomustine and vincristine (PCV), see e.g. [1,15,32,33] for the mechanism of action of those drugs. However, due to PCV's higher toxicity profile [7,24,29], TMZ is currently a preferred chemotherapeutic agent for LGGs' patients. Most commonly used TMZ treatment consists of doses of 150-200 mg per m 2 of patient body surface area administered in a sequence, or a "cycle", using the oncological nomenclature. The usual TMZ cycle takes 28 days and consists of doses given during the first five days of the cycle once per day, i.e. five days of drug administration are followed by 23 days of drug-free interval before the next cycle. However, such a treatment protocol has been based on a principle of "maximum tolerated dose over a short period of time" applied originally for high grade and fast-evolving tumours. Thus, new treatment regimens for LGGs are being investigated [17,22,48,49] in order to find the most effective chemotherapy schedule.
At the same time, there is a trend to administer chemotherapy using metronomic fractionation, i.e. schedules consisting of many, equally spaced and generally low doses of chemotherapeutic drugs without extended rest periods (see e.g. [2,3]). Thus, there are numerous, clinically relevant open questions regarding, among others, the proper timing and dosing of chemotherapy. These questions could be addressed using mathematical models constructed on the basis of biological and clinical observations and validated using patient data. Unfortunately, up to now, there have been very few mathematical models considering chemotherapy for LGGs. An interesting mathematical model addressing response of LGGs to TMZ was formulated by Ribba et al. in [43]. However, that model, due to its complexity, includes ten parameters and some of them are very difficult (or impossible) to measure or estimate based on available data, cf. [10]. Moreover, the model proposed in [43] was intended to reflect LGGs response not only to chemotherapy with TMZ, but also to chemotherapy with PCV and radiotherapy. However, all that phenomena was modelled in the same way. This seems to be a small oversimplification, as radiotherapy is well-known to induce a different response to chemotherapy and moreover there are clinically-recognised mathematical models capable of describing tumour response to radiotherapy sufficiently well, see [14,37]. The model of Ribba et al. was later extended by Mazzocco et al. in [30], where the authors took into account the possibility of acquired resistance to drug. That model was also used to find more efficient schemes for LGG chemotherapy with PCV in [31], however the authors drew their conclusions based solely on the results of simulations of that model.
In [5] we proposed a simple mathematical model for LGGs growth and response to TMZ. In order to describe the tumour growth and chemotherapy treatment we have introduced only three patient-specific parameters postulating that other parameters can be estimated from available studies for patients with LGGs and their response to treatment with TMZ [12,40,44]. The formulated model fits very well to longitudinal volumetric patient data from Bern University Hospital, for details see [5].
In this paper, we first describe a model proposed in [5], where a pharmacokinetic dynamics of the drug is modelled directly by an impulsive ordinary differential equation. Next, we generalise it by: (i) considering generalised tumour growth function and (ii) introducing into the system a new parameter γ > 0, which represents the effectiveness of the competition between the damaged cells and the proliferating ones.
We would like to stress that in [5] tumour growth was modelled by logistic function. However, we take into account a fact that the specific form of LGG growth function remains unknown. Thus, based on quantitative data from clinical studies (e.g. [35]), we only assume its general shape and in consequence we consider a broader set of functions that could describe the LGG evolution. We study the dynamics of resulting general model, showing that this dynamics is in some sense similar regardless of the choice of the specific growth function. Moreover, in [5] it was assumed that competition strengths of proliferating and damaged cells are equal. Biological intuition suggests that these strengths may not in fact be equal, hence we introduce the parameter γ described above. In [5] the function (describing the concentration of the drug) was assumed to be a solution to an impulsive equation with a finite number of impulses equal to the number of doses administered to a specific patient. In presented study we analytically investigate the effect of chemotherapy administered in cycles. However, in addition we also consider asymptotically periodic drug concentration function, which is a generalisation of the previously used one.
We would like to underline that results presented in [5] do not include mathematical analysis of the proposed model. Here, we investigate the mathematical properties of the proposed generalised model, including the existence and uniqueness of solutions, as well as the existence of the invariant set. We discuss the stability of existing steady states. We also indicate the condition under which the trivial steady state is asymptotically stable in the case of asymptotically periodic treatment. We also show that in some cases of periodic treatment there exist periodic solutions. We also provide estimations of suggested minimal doses for patients with LGGs based on data from Bern University Hospital described in [5].
1.1. Model formulation. The model proposed in [5] assumes that after the onset of chemotherapy treatment a LGG is composed of two type of cells: proliferative P and irreversibly damaged by chemotherapeutic drug D. Such a modelling approach was chosen in order to represent one of the most relevant and distinguishing feature of LGGs' response to chemotherapy, namely the observation that LGGs show a prolonged response to chemotherapy, sometimes lasting even years after it has ended [8,38,44]. The tumour growth in [5] was modelled through a logistic function with coefficient ρ, where its inverse gives an estimate of an average cell doubling time. Based on biological reports, such as [28], we assumed that irreversibly damaged tumour cells D try to enter mitotic cycle with the same probability as those active, but they die after a mean value of k such attempts. The mean number of cells damaged by the drug in a time unit was considered to be proportional to the concentration of the drug C multiplied by the average number of proliferating tumour cells with the rate α, describing the effect of TMZ acting on tumour cells. The model by Bogdańska et al. was formulated in a form of the following two ordinary differential equations, [5]: where P and D denote the mean number of LGG proliferating and irreversibly damaged tumour cells, respectively. In (1.1) function C describes the concentration of the chemotherapeutic agent (i.e. TMZ) and s ∈ R + 0 = [0, +∞). Clearly, because it is assumed that before s = 0 no treatment was administered, we restrict ourselves to the initial condition We considered a pharmacokinetic dynamics of the drug modelled directly by an impulsive ordinary differential equatioṅ where s j ≥ 0, j = 1, . . . , n, denote times of drug administration and C j are fractions of the dose administered at time s j that reach tumour tissue. Clearly, all model parameters are assumed to be positive due to their biological interpretation provided at the beginning of this Subsection. For finite number of chemotherapy doses the model dynamics is very simple. After the last dose administered at time s n the concentration of TMZ, C, decays exponentially to 0, thus all solutions to system (1.1) converge to the point (K, 0), as it was described in [5].
In this paper we propose a generalisation of the model described above, assuming that the tumour growth term is governed by a function having certain properties. Moreover, we consider that the chemotherapy drug is either administered in infinitely repeating cycles, or its concentration is constant. Thus, in general, we assume that the function C is non-negative and continuous, and we focus on the two particular cases: the function that is either asymptotically periodic (we give a precise meaning to this term later) or constant. Rescaling the variables in system (1.1) by taking with initial conditions: and where z is a non-negative, continuous function describing rescaled concentration of the chemotherapeutic drug. In this paper we consider a generalised version of system (1.5), that iṡ The function f describes the tumour growth, so we impose the following assumptions: = 0, and f is defined at zero by f (0) = 0, resulting in logistic-type (A3a) or Gompertz-type (A3b) growth.
Here, the new parameter γ > 0 represents the effectiveness of the damaged cells in competition with proliferating ones. Although we believe that assumption γ ≤ 1 would be biologically justified, we also consider γ > 1.
Sometimes it is convenient to consider the variable u and then system (1.7) takes the forṁ Mathematical analysis of the model.

2.1.
Existence and uniqueness of solutions, and existence of an invariant set. In this section we formulate and prove theorems guaranteeing existence and uniqueness of solutions to system (1.7) with any non-negative/positive initial conditions. Moreover, we show that if γ ≤ 1, then there exists a compact set invariant with respect to system (1.7).
Proof. The local existence of solutions to system (1.7) for every initial data from (R + 0 ) 2 follows from the continuity of the function of the right-hand side of system (1.7) in (R + 0 ) 2 . Moreover, the right-hand side of system (1.7) is a C 1 function in (R + 0 ) 2 \ {(0, 0)} or, if assumption (A3a) holds, even in a whole (R + 0 ) 2 . This fact implies the uniqueness of solutions in the respective set.
Note that in the case of function satisfying condition (A3b) we may loose uniqueness of solutions if x(0) = 0 and y(0) = 0 (e.g. if we take f (v) = 1/ √ v and z = 0).
Proof. Let x 0 = x(0) > 0 and y 0 = y(0) > 0. If assumption (A3a) holds, the non-negativity of solutions to system (1.7) follows from the form of the right-hand side of system (1.7). On the other hand, if (A3b) is satisfied, then for x sufficiently small we haveẋ > 0, thus x(t) > 0 and therefore y(t) > 0.
To prove the second part of the assertion, we consider variable u defined in (1.8) that fulfils (1.9b). Clearly, (1.9b) is a non-autonomous differential equation with x being a non-negative function of time, and it has unique solutions. Moreover, u ≡ f −1 (0) is a solution of (1.9b) for γ = 1. On the other hand, for γ < 1 we havė u| u=1 = (γ − 1)zx < 0. Thus, u(t) ≤ 1 for all t ≥ 0, which completes the proof.

2.2.
Existence and stability of the steady states in the case of constant z. First, consider the case of z being a constant function. With some abuse of the notation we write z(t) ≡ z. An easy calculation shows that system (1.7) has up to three non-negative steady states: Clearly, the positive steady state P 3 exists if z < f (0) (i.e. if either z < 1 or assumption (A3b) holds). Let us define 3. Let f satisfy assumptions (A1)-(A3) and z(t) ≡ z. Then the steady state Proof. Assume first, that assumption (A3b) holds. Then for x + γy < f −1 (z) we haveẋ > 0. Thus the solution to system (1.7) leaves the set (x, y) : x + γy < f −1 (z) which implies that P 1 and P 2 are unstable. On the other hand, if assumption (A3a) holds, the Jacobi matrix at P 1 and P 2 reads , respectively. Thus, the assertion (ii) follows as well as (i) for z = 1. Now, we prove that P 1 is locally stable if z = 1. Because z = 1, f is decreasing and f (0) = 1 (see (A2)-(A3a)). Consequently we haveẋ(t) ≤ 0 for all t ≥ 0. The function x is non-increasing and bounded from below, thus it has a limitx as t → +∞. We havex = 0 orx + γy = 1. However, the latter case is impossible if is negative for sufficiently large y. Thus P 1 is locally stable. Now, we study local stability of P 3 . The Jacobi matrix at P 3 reads where we have used the from of y * , the definition of δ (see (2.1)) and the fact that f is a decreasing function. Clearly, the stability of steady state P 3 depends on the value of trace and determinant of the Jacobi matrix, that is Consequently, the characteristic polynomial of the matrix J(P 3 ) equals Because κδz > 0, the steady state P 3 is asymptotically stable for z > (γ − 1) δ κ+γ . This inequality is equivalent to a condition κz + δ > γ(δ − z), which is fulfilled for either z ≥ δ or for γ < κz+δ δ−z . This completes the proof. We now determine conditions under which the type of steady state P 3 changes. Steady state P 3 can be either a node or a focus depending on the relation between the model parameters. Let us denote Theorem 2.4. Let f satisfy assumptions (A1)-(A3), z(t) ≡ z, and the steady state P 3 of system (1.7) exist. Then the steady state P 3 is a node in the following cases: (i) z > δ and (κ > κ 1 or 0 < γ < γ 1 ,) (ii) z < δ, • 0 < γ < γ 1 or (γ > γ 2 and κ 1 < κ).
Proof. We check when roots of (2.2) are real. To this end, we calculate the determinant ∆(γ) of (2.2) and we write it as a function of γ obtaining It is easy to see that ∆ is a quadratic polynomial of γ with two real roots γ 1 and γ 2 given by (2.4) and (2.5), respectively. Note also that κ > κ 1 implies that the coefficient of γ 2 is positive. Thus, a simple analysis of quadratic polynomial yields theorem assertions.
In Fig. 1, we graphically illustrate the assertion of Theorem 2.4 presenting the class and stability of the steady state P 3 depending on parameters κ and γ. Note that for κ = δ/z we have γ 1 = 0, while for κ = 4δ/z equality γ 1 = 1 holds. In addition, condition γ ≤ 1 yields stability of steady state P 3 (provided it exists).
Additionally, in Fig. 2 we present exemplary phase portraits for γ = 1 and f (x+y) = 1−x−y. In that case steady state P 1 is either a stable node (for z > 1) or Figure 1. Sketch presenting the class of steady state P 3 depending on considered cases: (left) z ≥ δ; (right) z < δ. Blue and white areas represent sets of parameters for which P 3 is node and focus, respectively. Dots denote region where P 3 is stable, no patternregion where P 3 is unstable. a saddle (z < 1). The steady state P 2 is a saddle independently of model parameters, while P 3 is locally stable as γ = 1 (compare Theorem 2.3). To be specific, P 3 is a stable node for κ > 4(1−z) z and a stable focus for κ < 4(1−z) z . In Figs. 3 and 4 we present exemplary phase portraits with nullclines for f (x + γy) = 1 − x − γy in a case when z < 1 and z > 1, respectively. Proof. To prove that system (1.7) has no periodic solutions we use Bendixson-Dulac criterion. We consider a variable u defined by (1.8) and we study system (1.9). Let us denote by F the right-hand side of system (1.9) and consider the function ψ(x, u) = 1 xf (u) as Dulac function. We calculate For γ ≤ 1 the sign of divergence div(ψF ) is negative almost everywhere in Ω. As a consequence, the system (1.9) does not have closed orbits in Ω.
Proof. The assertion of the theorem follows from Lemma 2.5, the Poincare-Bendixson theorem, together with local stability of the steady states.
2.3. The case of asymptotically periodic z. First, following [4], let us recall the definition.
Definition 2.7. We say that function z : R + 0 → R is an asymptotically periodic function with a period T if there exist functions z per , z rest : R + 0 → R such that z per is periodic with a period T , lim t→+∞ z rest (t) = 0 and z(t) = z per (t) + z rest (t).
We show later that solution to (1.3) is an asymptotically periodic function with period T (in the sense of the above definition) if there exists p ∈ N such that s j+p = s j + T for all j ∈ N. Now, we show that for any asymptotically periodic function z whose mean value converges to a value not smaller than 1, the solutions to system (1.7) tend to the steady state P 1 . converges to P 1 .
Proof. Due to the fact that f is decreasing, γ > 0, x and y are non-negative we have and as a consequencė holds. The Gronwall inequality implies (2.6) We divide the integrand into three components. The first one includes function depending only on x, the second -periodic function defined on the interval whose length is a multiplication of the period, and in the last component we group the remaining terms. To this end, for any t ≥ 0 we find m ∈ N ∪ {0} andt ∈ [0, T ) such that t =t + mT and we rewrite (2.6) in the following manner where

MODEL OF CHEMOTHERAPY FOR LGG 2159
Note that R(t) is uniformly bounded becauset ∈ [0, T ) and the integral +∞ 0 z rest (s) ds is convergent.
First, considerz > 1. Then A(t) ≤ 1 as f is decreasing and x(t) ≥ 0 for all times t. Moreover, Thus, x(t) → 0 as t → +∞. Now, assume thatz = 1. Then B m = 1. As x(t) ≥ 0 for all t ≥ 0 and f is As a consequence, x(t) → 0 as t → +∞ due to the continuity of f . On the other hand, the divergence of the integral implies A(t) → 0 as t → +∞. This proves that x(t) → 0 as t → +∞ because B m is uniformly bounded, as stated before. This completes the proof.
The model behaviour forz < 1 is more complex. Numerical simulations suggest that the steady state P 1 is repulsive in that case, see an example in Fig. 5. However, if z is a periodic function and its values remain inside the interval (0, 1), then we can prove the existence of a periodic solution.  Proof. For convenience we consider system (1.9), that is system (1.7) in variables x and u. Moreover, Theorem 2.2 allows us to restrict our attention to the set For any (x 0 , u 0 ) ∈ Ω u we define the operator where x(t; x 0 , u 0 ), u(t; x 0 , u 0 ) is a solution to system (1.9) with initial condition ∈ Ω u is a fixed point for the operator Φ T : Ω u → Ω u , then, due to the assumption of periodicity of z, a solution to system (1.9) with initial condition (x 0 , u 0 ) is periodic.
We use the Brouwer theorem [6] to show that such fixed point of the operator Φ T exists. To this end, we show that there exists a convex compact subset K ⊂ Ω u such that Φ T (K) ⊂ K.
Vector field corresponding to the system (1.9) is the following one: Let us divide the set Ω u in four sets Now, we define the border of the set K. Let ϕ 1 be a trajectory of system (1.9) with z ≡ z m and arbitrary chosen initial condition (x 0 , x 0 ) ∈ Ω 1 for time [0, t * ], where t * is the first point t > 0 at which this trajectory reaches the boundary of Ω 1 . Because of the direction of the vector field F zm (which is independent of time) in Ω 1 , the trajectory ϕ 1 reaches the line (κ + γ)x = κu at t * . Let us call this point (x 1 , u 1 ). Let ϕ 2 be a trajectory of system (1.9) with z ≡ z M and initial condition (x 1 , u 1 ) for time [0, t * * ], where t * * is the first point t > 0 such that the trajectory ϕ 2 reaches the boundary of Ω 2 . Because of the direction of the vector field F z M (which is independent of time) in Ω 2 , the trajectory ϕ 2 reaches the line u = 1 − z M . Denote this point (x 2 , u 2 ). As ϕ 3 we denote the line that connects points (x 2 , u 2 ) with (x 2 , u 3 ) lying on the line κu = (κ + γ)x, thus u 3 = (1 + γk)x 2 . We denote by ϕ 4 the horizontal line between (x 2 , u 3 ) and (u 3 , u 3 ). Finally, we denote the line connecting points (u 3 , u 3 ) and (x 0 , x 0 ) by ϕ 5 . Obviously, the curve Γ = ϕ 1 ∪ ϕ 2 ∪ ϕ 3 ∪ ϕ 4 ∪ ϕ 5 appoints a convex set in Ω u . We denote the set bounded by this curve (with the boundary) by K, compare Fig. 6. We show, that if (x(0), u(0)) ∈ K then for all t ≥ 0 (x(t), u(t)) ∈ K, where (x(t), u(t)) is a solution to system (1.9). Two parts of the boundary of K are the trajectories of system (1.9) with z ≡ z * , where z * = z m or z * = z M . The normal vector to such trajectory pointing towards the interior of K is Let (x(t), u(t)) be a solution to system (1.9) and assume that for some t * it reaches the set ϕ 1 ∪ ϕ 2 . Then, the scalar product of the vector tangent to the trajectory with the normal vector n is Assume now, that (x(t * ), u(t * )) ∈ ϕ 1 . Then, z * = z m , thus z(t) ≥ z m and (κ + γ)x − κu ≤ 0 as ϕ 1 ⊂ Ω 1 . Thus, F z (t), n ≥ 0. Similarly, if (x(t * ), u(t * )) ∈ ϕ 2 then z * = z M , hence z(t) ≤ z M and (κ + γ)x − κu ≥ 0 as ϕ 2 ⊂ Ω 2 . Thus, again F z (t), n ≥ 0. Note, that in Ω 3 we haveẋ(t) > 0. Similarly, in Ω 4 the derivative of u is positive. Thus, on the lines ϕ 3 and ϕ 4 , the vectors of the vector field points towards the interior of K. Finally, it is easy to see that on ϕ 5 the vector field also points towards the interior of K because, as it was shown before, solutions x(t), y(t) are positive and u(t) = x(t) + γy(t) ≤ 1. This implies that solution to system (1.9) that starts in K cannot leave this set. And therefore, Φ T (K) ⊂ K. In Fig. 6 we present the sketch of the set K, curves ϕ i , i = 1, . . . , 5, and sets Ω j , j = 1, 2, 3, 4.
Thus, the Brouwer theorem yields existence of a fixed point of Φ T and, in consequence, a periodic solution to system (1.7).
2.4. Chemotherapy given in cycles. As mentioned before, variable z represents the rescaled concentration of the chemotherapy drug C. The dependence between those two variables is given by (1.4). In this section, we relate theoretical results from the previous subsections with medical data. To this end, we use non-scaled form of the function C describing the concentration of the drug. If C andC are solutions to (1.3) with jumps C j andC j at moments s j ands j , j ∈ N, respectively, then the function C +C is a solution to (1.3) with jumps C j , C j at moments s j ,s j , j ∈ N. Proof. It is an easy adaptation of the result from Section 3.1 in [39]. Proposition 2.11. Assume that the impulse times s j and jumps C j , j ∈ N, are such that there exist T > 0 and p ∈ N such that s j+p = s j + T, and C j+p = C j .
Then the solution to (1.3) with such periodic impulses is an asymptotically periodic function C(s) = C per (s) + C rest (s) with Proof. It is easy to see that if impulse times are s j = s 0 + jT , then assertion of Proposition 2.10 does not change except for the formula for C rest (s), which now reads Because the impulses are periodic, we divide impulse moments into p subsequences (r j n ) n∈N , j = 1, . . . , p such that r j n = s j + (n − 1)T, compare with Fig. 7. Moreover, for each sequence (r j n ) n∈N , the corresponding jumps are the same and equal to C j . Combining this with assertion of Proposition 2.10 and Remark 1, we deduce the assertion of Proposition 2.11. Corollary 2.12. If the assumptions of Proposition 2.11 hold, then a mean valueC of the function C per does not depend on the impulse moments, but only on the total sum of jumps and the length of the period. This mean value is given by the formulā C = C 1 + C 2 + · · · + C p λT .

MODEL OF CHEMOTHERAPY FOR LGG 2163
The above results mean that the asymptotic mean value of the treatment depends only on the clearance rate λ and the mean drug administration per time unit. In particular, it does not depend on the distribution of doses, nor on the dose timing.
3. Discussion. In the presented paper we study a model of LGG growth and its response to chemotherapy. The considered model is a generalisation of the model proposed in [5], where the particular (logistic) form of the tumour growth function was considered. Here, we take into account the possibility that a more general tumour growth function may fit better to the individual patient data. Thus, we only assume its general properties (determining its shape) and consequently consider a broader class of possible functions. In addition, we incorporate a new component to tumour growth term -a parameter γ describing the effectiveness of the competition between damaged and proliferating cells. We prove the existence and uniqueness of solutions together with the existence of the invariant set. We study the long time behaviour of the model depending on its parameters, including the ones describing chemotherapy effect. We discus the stability of the steady states and we prove the condition under which the trivial steady state is asymptotically stable in the case of periodic treatment. We also show that there exist periodic solutions for some specific cases of periodic treatment function.
Note that in the chemotherapy scheme most frequently used for LGG patients, equal doses of TMZ (namely 150 mg per m 2 of body surface area) are administered. However, it has been noted both in clinical and mathematical studies [17,22,34,43,45,48] that such a dose seems not to be the most appropriate for slowly-growing LGGs. Our generalised model given by system (1.7) together with Theorem 2.8 yields a few possible therapeutical implications. It suggests the possibility of estimating the minimal dose d allowing to eradicate tumour under assumption that drug doses are administered infinitely many times and that tumour cells do not acquire drug resistance. Clearly, in real life one should focus on treatment schemes with finite time horizons. However, a significantly prolonged chemotherapies could be considered such as therapies lasting until obtaining the expected response to treatment. In [18], Khasraw et al. report cases of LGG patients for whom chemotherapy treatment administered for even 5-8 years did not cause serious side effects. Thus, on the basis of this clinical study, we presume that prolonged chemotherapy could remain safe for LGG patients. In addition, Mannas et al. conclude that longterm chemotherapy with TMZ could be considered a therapeutic option as long as appropriate monitoring is assured [27]. Now, we use theory developed in the paper to estimate minimal effective dose in the case of long time treatments. In order to do that, we recall the values of model parameters. In [5] the parameter λ was fixed to be 8.3184 per day, while parameters α and ρ describing the effectiveness of the drug and the proliferation rates of LGG, respectively, were assumed to be patient-specific and were fitted for each patient independently. In [5] we also estimated that a single fraction of drug dose actually acting on tumour tissue was proportional to the administered dose with coefficient ν=4·10 −6 m 2 /ml. Finally, considering model rescaling (1.4) we havē patient id minimal drug dose (in mg) Figure 8. The minimal eradication dose d per m 2 of body surface estimated for the patients' parameters estimated in [5]. estimated for patients considered in [5], see also Fig. 8. Clearly, dose of 5.42 mg is much smaller than the dose of 150-200 mg that is administered in a clinical practice. This result suggests that so called metronomic therapy could be a good option for successful treatment of selected patients.
There have already been several studies analysing metronomic chemotherapy for HGGs, see e.g. [21,49]. As for LGGs, Lashkari et al. in [22] studied two alternative treatments with smaller TMZ dose. The authors concluded that metronomic regimens of TMZ may result in better LGGs response compared to the conventional regimen, however there is a need for randomised clinical trials to verify this result. Also, recent observations of the benefits of metronomic chemotherapy with other drug (vinblastine) administered for a limited number of LGG patients may be related to our analysis [13]. Our result also indicates that more patient specific therapy schemes should be designed to obtain the best effectiveness. It also points that the proper dosing should be rather decided on the basis of tumour-specific characteristics.
We would also like to emphasise that doses smaller than those calculated from (3.1) could not be effective in treating gliomas in an infinite time horizon, thus they should not be considered for realistic finite time schemes.
On the other hand, in reality one also needs to take into account the possibility of drug resistance that can be acquired by tumour which lies beyond presented considerations. Thus, in the future, we plan to extend the model in order to be able to broaden the possible applicability of this results and address other clinicallydriven problems such as drug resistance and cytoxicity being a strongly limiting factor to any chemotherapy. Such a task would require inclusion of more terms and variables in the model and more biomedical data in order to validate the model and estimate model' parameters.
In particular, the main cytotoxic effects induced by TMZ are due to haematologic toxicity [1,41,46,47]. In order to include them in considered model, one would require additional data from systematic blood morphological patients tests, preferably before and during the treatment with TMZ. Moreover, in [19,20] it has been discussed that TMZ effectiveness in treating gliomas is limited not only due to cytotoxicity, but also due to the drug influence on endothelial cells. If we wish to include that phenomenon in our model as well, we would need to consider additional variable in the system. Consequently, the model should be calibrated with additional data. Unfortunately, at the moment we do not have access to patients TMZ cytotoxicity data, nor the quantitative data on this specific drug effect on endothelial cells.
The effect of acquiring resistance can be modelled in various ways, see e.g. [9,11,23,36] and references therein. However, the selection of suitable approach needs to be deeply investigated, mainly by taking into account available data and the possibility of verification of the obtained results. To conclude, we plan further studies aimed to find the best possible way to build and validate models describing the above phenomena in more detail which would also be used to address questions concerning LGG which remain open.