Dynamical Systems Properties of a Mathematical Model for the Treatment of CML

A mathematical model for the treatment of chronic myeloid leukemia (CML) through a combination of tyrosine kinase inhibitors and immunomodulatory therapies is analyzed as a dynamical system for the case of constant drug concentrations. Equilibria and their stability are determined and it is shown that, depending on the parameter values, the model exhibits a variety of behaviors which resemble the chronic, accelerated and blast phases typical of the disease. This work provides qualitative insights into the system which should be useful for understanding the interaction between CML and the therapies considered here.


Introduction
Chronic myeloid leukemia (CML) is a hematologic cancer that accounts for about 15% of all leukemias in adults and is characterized by uncontrolled expansion of myeloid cells in the bone marrow and their accumulation in the blood [1].The progression of the disease can be divided into three phases denoted chronic, accelerated and blast [1].The chronic phase can last several years with levels of immature white blood cells (blasts) growing steadily but at a low rate.Once the disease enters the accelerated or blast phase, cells proliferate rapidly and the disease can be lethal within a few months if not treated.Current standard of care includes targeted tyrosine kinase inhibitors (TKIs), which have significantly improved long-term survival rates [2].
Responses to certain treatments have offered evidence of an immune component in the disease [3].Early indications were provided by a correlation between incidence of graft-vs-host disease and improved leukemia-free survival in CML patients who had received allogeneic stem cell transplants [4].Additionally, treatment with interferons (which are known to be immunomodulatory) has led to complete or partial responses in some fraction of CML patients [5].More recently, studies that include immunomodulatory therapies such as nivolumab have been initiated [6].
Mathematical modeling of CML dynamics has a history dating back to the late 1960s with early work of Rubinow and Lebowitz [7,8].Models by Fokas et al. [9] in the 1990s focused on maturation and proliferation of T-cell precursors.In 2004, Moore and Li [10] published a model of CML dynamics, which accounts for the actions of naive and effector T-cells separately.In [11], this model was analyzed as an optimal control problem.The model presented here first appeared in [12] and models the immune system effects with one compartment, and separates the CML cells into quiescent and proliferating classes.The rationale behind this new model is the ability to represent certain types of therapies for use in combination treatment.These therapies are: a BCR-ABL1 tyrosine kinase inhibitor (e.g., a therapy such as imatinib), an immunomodulatory therapy (e.g., a therapy such as nivolumab), and a therapy that combines both actions (e.g., a therapy such as dasatinib).
The model introduced in [12] is reviewed in Section 2 and then analyzed as a dynamical system with constant drug concentrations in Section 3. The analysis is carried out theoretically for values of parameters covering a range of dynamic possibilities.As will be seen, there are parameter values for which the model can have an asymptotically stable equilibrium point in which all the state variables are positive.This could be interpreted as disease control through continuous therapy.As parameters change, the system can become unstable and undergo exponential growth, representing the accelerated or blast phases of the disease.Our analysis incorporates constant drug concentrations, and thus provide insights into the dynamics both without and with treatment.In particular, we analyze how an increase in the levels of each of the three treatments affects the values of all three populations, the two types of leukemia cells and the strength of the immune effect.The combination of theoretical analysis and simulations is intended to shed some light on understanding the long-term dynamics of this disease under treatment.

A Mathematical Model for the Treatment of CML with BCR-ABL1-Targeted and Immunomodulatory Drugs
The mathematical model below was originally published in [12] in 2015.

A Brief Review of the Mathematical Model
Let Q be the concentration of quiescent leukemic cells, P the concentration of proliferating leukemic cells, and E the strength of immune system effects.We will consider E to represent effector T cell concentration levels, and will refer to E in the remainder as a concentration of effector T cells.The model contains three controls u 1 , u 2 and u 3 that all denote normalized levels of different therapies.The roles of the specific drugs are illustrated in Figure 1 taken from [12] with arrows indicating amplification of effects and vertical bars indicating inhibition.The control u 1 represents the normalized concentration of a BCR-ABL1 inhibitor (such as imatinib) that mainly has an inhibitory effect on the highly-proliferating leukemic cells; u 2 is a BCR-ABL1 inhibitor that inhibits BCR-ABL1 that also has immune effects (such as dasatinib); while u 3 represents an immunomodulatory compound (such as nivolumab).In this system, all parameters are non-negative.For P = 0 or E = 0, we extend the system by defining it using the limits as P → 0 or E → 0, respectively.The cell count numbers for Q are relatively small and are therefore modeled by an exponential function with growth coefficient r Q .For the proliferating cells P we model growth with a Gompertz function, as Afenya and Calderón state that this is best for describing CML growth [13].The immune effect E (effector T cells) also has its rate of increase modeled by a Gompertz function, so as to have approximately exponential growth when numbers are very small, but still be bounded above.In the populations P and E, replication rate constants are represented by r P and s E , and carrying capacities (or steady states) by P ss and E ss , respectively.The natural death rate constants of the respective populations are denoted by δ Q , δ P and δ E .The population Q consists of leukemic cells that are quiescent.Some or all of quiescent leukemic cells may be stem cells [14].When quiescent cells divide, one copy is assumed to be the same kind as the original cell while the second copy may differentiate further into a proliferating type.For this reason, the transition term k P Q is not subtracted from the quiescent cell population in (1).This term represents the rate at which quiescent cells produce differentiated proliferating cancer cells, with the population Q the source for the population P.
The control variables represent the concentrations of the respective drugs, and their effects (pharmacodynamics) are modeled by Michaelis-Menten terms with different maximum effectiveness on the various populations.In modeling the combined drug actions it is assumed that any two drugs act independently of each other.Thus the term represents the effects that drugs 1 and 2 have on decreasing the proliferation of the population P.
A term of the type represents the enhancement of the actions of the effector T cells E on the quiescent cells Q as a consequence of the activities of drugs 2 and 3.In each of the equations, the enhancement and inhibition effects of the drugs by means of the immune system are modeled additively.
The "C 50 " parameters U1C 50 , U2C 50 , and U3C 50 represent the concentrations required to achieve half of the maximal effects of u 1 , u 2 , and u 3 , respectively.These and EC 50 and PC 50 are assumed to be fixed across effects being modeled.These represent "potency" levels depending intrinsically on the particular therapy or population, and not on the setting of the effect.The maximum possible effect size is allowed to depend on the setting.
The equations above represent a semi-mechanistic, fit-for-purpose, minimal model.It is minimal in the sense that it only includes the levels of cell interactions needed to allow the controls to have their expected effects.Some of the terms are based on models validated with data, but other terms take forms that are more heuristic.For example, all of the control effect terms take a Michaelis-Menten or "Emax" form.This is because we wish to model very small effect at low levels of drug, as well as a limiting or asymptotic maximal effect at high levels of drugs.We chose the simplest among the models with this behavior that are typically used in drug development [15].
The states, controls, and related parameters are listed in Tables 1 and 2. Table 1 gives those parameters that are unrelated to the drug actions and make up the untreated, or uncontrolled, system; Table 2 lists the treatment-specific parameters in the model.In this paper, we do not fit or fix specific parameter values, and instead analyze the dynamic properties of the system (1)-(3) for large ranges of possible values.We include in the tables below two different sets of numerical values that we use to illustrate the dynamic properties of the system.These parameter values are purely for numerical illustration and do not reflect specific model fits or therapies.The focus of this paper is the mathematical analysis of the entire system rather than an analysis for particular parameter values.maximum possible effect of u 3 on 5 0.05 stimulating proliferation of effector T cells maximum possible effect of u 3 on 0.7 0.07 prevention of the death of effector T cells U3C 50 concentration of u 3 that gives half the maximum effect 0.7 0.7

Scaling of Parameters
We note that the dynamical system has various groups of symmetries that can be used to scale the variables and controls.Here we normalize all the "C 50 " parameter values to 1 by rescaling the corresponding variables in terms of these quantities.This simply minimizes the number of parameters to be considered in the analysis of the system.For example, let Q ref be a constant to be determined later, and define and Then we have that and analogously for the other terms.
For the differential equations, we obtain Under this scaling all remaining parameters in this equation are invariant and need not be changed.Similarly, and Thus, if we re-scale and the steady-state values as then formally the equations are the same as before with all "C 50 " values in the Michaelis-Menten expressions normalized to 1.All other parameters remain unchanged and even their interpretation is the same as before.For the theoretical analysis and numerical computations this eliminates five parameters and introduces a favorable scaling to the variables.Naturally, the original parameters are still calculated for an interpretation of the results.

System Properties for Constant Concentrations
CML has three distinct phases, a chronic one that can last from three to five years, during which leukemic cell counts are low but may grow steadily, and accelerated and blast phases that may last for a only a few months and are characterized by higher cell counts or a rapid increase in cell counts followed by death of the patient [10].Here we analyze the dynamical system to determine if it can capture such features.

Reduction to the Uncontrolled System and Basic Dynamical System Properties
We carry out the dynamical systems analysis for constant controls, i.e., concentrations.We do not explicitly include pharmacokinetics (fluctuations in concentrations that depend on doses).The treatments considered are either administered daily or have long half-lives, and such pharmacokinetics are not expected to be significant for the treatment periods we consider here (five years or longer).We also mention the 2009 paper by Shudo et al. [16] that supports this assumption in the setting of hepatitis C.
Keeping the "C 50 " parameters in their original formulation in the controls, we define new drug-dependent parameters as With these identifications, the dynamical system with constant controls is identical with the uncontrolled system and therefore, without loss of generality, the analysis can be done on the uncontrolled system.Returning to the original notation without the carets, we thus consider the following equations: The model with an exponential growth term on Q has various long-term behaviors.These include the extremes in which Q decays exponentially to zero or grows exponentially beyond limits, but there also is the possibility that nontrivial equilibrium points (Q * , P * , E * ) exist for which all three populations are positive.The first case corresponds to a scenario in which the patient goes into a stable deep molecular response.For the uncontrolled system, this may not seem to be of interest, but since the model includes the case with controls, this gives us information about which combinations of constant concentrations of the drugs would lead to an eradication of Q.The case of exponential growth may characterize the accelerated or blast phase as these phases have short doubling times [17].The conditions under which this is the long-term behavior of the system give information about what controls are needed for successful treatment.An asymptotically stable equilibrium point (Q * , P * , E * ) with positive values could be interpreted as describing a subset of the chronic phase where net growth rate is zero, controlled by therapy or immune effects.Depending on the values of the parameters, this equilibrium point may be stable or unstable.Since in real life parameters may not be constant, bifurcation phenomena would be a mathematical description of the transition from chronic to the accelerated or blast phases.Knowing the parameter values when this may occur would be of interest.Our aim in the following is thus to determine the asymptotic behavior of the trajectories of the system.
We start with some basic properties.The positive orthant is positively invariant for the dynamics.This is because the planes Q = 0 and E = 0 are invariant under Equations ( 6) and ( 8) and Ṗ ≥ k P Q whenever P = 0. Thus, starting at a positive initial condition (Q 0 , P 0 , E 0 ), it follows that the solutions remain positive for all times.For the long-term behavior of the system, the equilibrium solutions in the closure of P, P = clos(P), also matter.Recall that the system is defined and continuous on P due to the use of the limits as P → 0 and E → 0 in place of P = 0 and E = 0, respectively.The vector field defining the P and E dynamics is not continuously differentiable at P = 0 or E = 0, but these values are repelling and thus this does not become an issue.
Lemma 1.The equilibrium solution E * = 0 is repelling: there exists a positive threshold E ∆ < E ss such that dE dt is positive on (0, E ∆ ].In particular, once E ss > E(τ) ≥ E ∆ , then E(t) ≥ E ∆ for all t ≥ τ.Furthermore, for E(0) < E ss , E will remain below E ss .
Proof.The terms in the last parentheses in Equation ( 8) are bounded between 1 and 1 + P max,2 and thus, as E → 0, the Gompertzian growth dominates the dynamics.Specifically, let Lemma 2. The equilibrium solution P * ≡ 0 is repelling: there exists a positive threshold P ∆ < P ss such that dP dt is positive on (0, P ∆ ].In particular, once P ss > P(τ) ≥ P ∆ , we have P(t) ≥ P ∆ for all t ≥ τ.
Proof.For values of E less than E ss , we have that This proves the result.
Note, however, that P is not necessarily bounded.For, with P = P ss , Equation (7) becomes and thus, if Q is large enough, this term will be positive.Hence, if Q grows exponentially, P will diverge to +∞.
Proof.We need to show that for every positive value P there exists a time T so that P(t) ≥ P for all t ≥ T. We first remark that P is unbounded.For, if there exists a value P with P ss < P < ∞ so that P(t) ≤ P for all times t, then the term r P ln P ss P is bounded below.
By assumption, there exist positive constants α and β so that Q(t) ≥ αe βt for all t.Hence, for t sufficiently large we have that dP dt (t) = r P ln P ss P(t) Contradiction.
Given P ≥ P ss , choose Ť so that Since P is not bounded, there exists a first time T > Ť so that P( T) = P + 1.We claim that P(t) > P for all t ≥ T. For, if there exists a time τ > T such that P(τ) = P, then dP dt (τ) = r P ln > r P ln Contradiction.Thus P diverges to +∞.
3.2.Dynamics on the Plane Q = 0 The plane Q = 0 is invariant under the dynamics and can have regions that are repelling or attractive.We first analyze the reduced dynamical system in this boundary stratum of P, i.e., consider the equations and denote by P0 its closure, P0 = {(P, E) : 0 ≤ P ≤ P ss , 0 ≤ E ≤ E ss }.For Q ≡ 0 the variable P is bounded above by P ss and therefore the compact set P0 is positively invariant under Equations ( 9) and ( 11).The dynamical system has the following trivial equilibrium solutions in the boundary of P0 : (0, 0), (P * , 0) with P * = P ss exp − δ P r P < P ss and (0, E * ) with E * given by In view of Lemmas 1 and 2 these solutions are unstable.While the origin has two unstable modes, the equilibrium points (P * , 0) and (0, E * ) are saddles with the respective axes forming the stable manifolds and the unstable modes entering the interior of P 0 .It is clear from this that there needs to exist at least one more equilibrium point (P * , E * ) in P 0 .Lemma 4.There are no periodic orbits in P 0 .Then equilibrium values P * are fixed points of the function Φ, P = Φ(P), in the interval [0, P ss ].Since Φ(0) > 0, Φ(P ss ) < P ss , and Φ is continuous in P, it follows that there exists at least one solution.The derivative Φ of Φ is given by

Proof
and thus has the same sign as Ξ (P).Now Ξ (P) = P max,2 − P max,1 (1 Thus Φ is strictly increasing for P max,2 > P max,1 and strictly decreasing for P max,2 < P max,1 .If P max,2 = P max,1 , then Equilibria are intersections of the graph of Φ with the diagonal and thus there exists a unique equilibrium point (P * , E * ) ∈ P 0 if P max,1 ≥ P max,2 , but multiple solutions are possible if P max,1 < P max,2 .
We determine the stability of (P * , E * ) for the reduced system, i.e., within the invariant plane Q = 0.The Jacobian matrix at the equilibrium point is given by Using the equilibrium relations we can write the (2, 1)-term as .
If we write χ(t) = t 2 + a 1 t + a 0 , then a 1 is positive and thus the equilibrium point is locally asymptotically stable if a 0 is positive while it is unstable if a 0 is negative.A saddle node bifurcation occurs as a 0 = 0.It immediately follows that (P * , E * ) is locally asymptotically stable if P max,1 ≥ P max,2 , i.e., if the stimulating effect of the tumor on the effector cells is larger than the inhibiting effect of the tumor on the effector cells.We have the following result: Proposition 1.If P max,1 ≥ P max,2 , then there exists a unique equilibrium point (P * , E * ) in P 0 and it is globally asymptotically stable in the sense that its region of attraction is the full rectangle P 0 .
Proof.The set P 0 is positively invariant and every trajectory γ starting in P 0 has a non-empty ω-limit set Ω(γ).Because of the stability properties of the equilibria in the boundary of P 0 , this ω-limit set Ω(γ) lies in P 0 .Since there exist no periodic orbits and since (P * , E * ) is the only equilibrium point, it follows from Poincaré-Bendixson theory that Ω(γ) = {(P * , E * )}, i.e., all trajectories starting in P 0 converge to (P * , E * ) as t → ∞.
It is clear from Poincaré-Bendixson theory that even if P max,1 < P max,2 , the equilibrium point (P * , E * ) is globally asymptotically stable (in the sense that its region of attraction contains the set P 0 , and only this region is relevant for the problem) as long as it is the only equilibrium point in P 0 .This is shown in the phase portraits for the uncontrolled system in Figure 2; Figure 3 shows a case where P max,1 > P max,2 .(The values of the parameters are given in Tables 1 and 2.) We also show the phase-portraits for the systems when one of the controls is set to be equal to 1 and all others are zero.The two sets of figures illustrate two different scenarios, one where the control parameters are such that the equilibrium can be effectively controlled by all the drugs (Figure 2), the other where it is essentially only the control u 1 that is able to move the equilibrium point.However, this behavior depends on the fact that Q = 0.  1 and 2. Phase portraits of the reduced dynamics for Q = 0 and P max,1 > P max,2 for the uncontrolled system (top, left) and for constant controls u 1 ≡ 1 (top, right), u 2 ≡ 1 (bottom, left) and u 3 ≡ 1 (bottom, right).The numerical values for these phase portraits are given in Tables 1 and 2.
Since the coefficient a 1 is always positive, as a 0 vanishes the Jacobian matrix has the eigenvalue 0 and the other eigenvalue is negative.At such a point saddle-node bifurcations arise and two new equilibria, one stable, the other unstable, are born. .
This condition is equivalent to (14).
For the underlying biological problem it is natural that an inhibition effect would be smaller than a stimulation effect.Also, the denominators are quadratic in the respective variables E and P, but these variables are scaled.In principle it is possible to satisfy ( 14), but we did not come across this in our simulations.

Dynamic Behavior for Positive Q-Values
For the behavior of the overall system, the Q dynamics are essential.If one considers the above equilibria in the plane Q = 0 now in the full three-dimensional space, then the first row of the Jacobian matrix at (0, P * , E * ) takes the form while the local stability properties for the overall system are the same as in the (P, , then there exists a 1-dimensional center manifold (corresponding to the 0 eigenvalue).In this case we have Q = 0 and there exists a curve of equilibria emerging from (0, P * , E * ) parameterized by Q or P (also see below).
Generally, (1)-( 3) is a time-varying linear system dominated by exponential growth and decay, depending on the parameter values.If then Q grows exponentially and no steady state exists.In this case, the influx k P Q eventually becomes the dominant term in Equation ( 2) and P also grows beyond limits (Lemma 3).This represents the malignant scenario in the model which corresponds to a highly-aggressive form of the disease or the accelerated or blast phase.The other extreme arises if r Q < δ Q .In this case Q exponentially decays to 0 for the uncontrolled system and overall trajectories converge to one of the equilibria (0, P * , E * ) in the plane Q = 0.If there exist multiple such equilibria, there exists a stable manifold for the unstable one that separates the regions of attraction for the stable equilibria.This would reflect a scenario when Q initiates the disease, but eventually dies off and the remaining P population determines the outcome of the disease.This could be benign if P * is small (a form of successful immune surveillance) or malignant if this value is larger.In such a case, however, one only needs to deal with the proliferating cells as far as treatment is concerned.This appears less likely (unless it could be induced by the drugs) and in the uncontrolled case of the disease we would have r Q > δ Q .
The interesting and most difficult case arises when the uncontrolled system has a chronic steady state or undergoes exponential growth without treatment, but has a negative net growth rate for Q with treatment.This is the case if the parameters satisfy the following condition (A): or, with the in the original form, Thus the replication rate constant r Q needs to be greater than the death rate constant δ Q (this naturally will be satisfied for parameters in a disease state), but at the same time, the drugs need to be able to raise the maximum effectiveness Êmax,1 = E max,1 1 U3C 50 +u 3 high enough that the magnitude of the immune system effect can overcome the difference.These appear to be natural conditions.Assuming that (15) holds, there exists a unique value E * ∈ (0, E ss ) for which Q = 0, namely with Q increasing for E < E * and decreasing for E > E * .In this case, the interplay between the variables allows for a steady state (Q * , P * , E * ) to exist with all values positive.We call such an equilibrium point (Q * , P * , E * ) positive.

Special Case: P
We first discuss the dynamical behavior of the system for the case P max,1 = P max,2 which is quite different from the cases P max,1 = P max,2 .If these effective rates are equal, we have that and it follows that E is strictly increasing for E < E * = E ss exp − δ E s E and strictly decreasing for E > E * .Therefore, as t → ∞, the E-dynamics approach E * , monotonically increasing if the initial condition is smaller, monotonically decreasing if it is higher.Consequently also the Q-dynamics approach the steady-state behavior and decrease exponentially if In the first case this also generates unbounded growth in P (Lemma 3) leading to behavior consistent with the blast phase of the system.In the second case, Q decays exponentially to 0 and P converges to the unique and asymptotically stable equilibrium point P * on Q = 0. Overall, and writing in the constant controls (the respective concentrations u i ) we have the following result: Proposition 3. Suppose P max,1 = P max,2 and let then all trajectories (Q(t), P(t), E(t)) converge to the unique and asymptotically stable equilibrium point (0, P, Ê) in the boundary of P 0 , whereas if then Q grows exponentially and lim t→∞ P(t) = +∞ and lim t→∞ E(t) = Ê.If We analyze whether positive equilibrium points (Q * , P * , E * ) exist.Throughout this section we assume that condition (15) is satisfied, i.e., that since otherwise Q grows exponentially.
Proof.The equilibrium relation for Equation (6) uniquely determines E * : Given E * , the equilibrium condition on the effector cells, Ė = 0, is equivalent to The quantity (17) only has the solution P * = 0; otherwise there exists a unique solution P * = P * (E * ) given by Using the equilibrium relation for E * , this can equivalently be expressed in the form Note that Q * is positive if P * ≥ P ss while otherwise this becomes a requirement on the equilibrium value P * = P * (E * ), namely If E max,2 = E max,1 , then this simply becomes P * > P ss exp − δ P r P r Q δ Q . In either case, there exists at most one positive equilibrium point given by Equations ( 16), ( 18) and (20).
Remark 1.As P max,1 → P max,2 , condition (15) implies that along a positive solution P * (E * ) we must have and thus the limit taken along these positive solutions only exists if E * → Ê = E ss exp − δ E s E and if In this degenerate case, the equilibrium conditions Q = 0 and Ė = 0 are automatically satisfied and there exists a one-dimensional equilibrium manifold, namely M = Q * (P), P, Ê : P > 0 with the P value arbitrary and Q * (P) given by We now investigate the stability of the positive equilibrium point.The partial derivatives of the equations defining the dynamics at the equilibrium point are given by Note that the equilibrium condition for The characteristic polynomial for the Jacobian matrix is given by .
If P max,1 < P max,2 , then a 0 is negative and the positive equilibrium point is unstable, i.e., once the maximal inhibiting effect of the tumor on the effector cells is larger than the maximal stimulating effect, no steady-state positive solution exists.Note further that for a 0 < 0 the characteristic polynomial χ(t) = t 3 + a 2 t 2 + a 1 t + a 0 has exactly one change of sign in its coefficients and thus there exists a unique positive root.So the equilibrium point has a two-dimensional stable manifold that separates the regions where Q and P diverge to infinity from the region where Q converges to 0. Thus we have the following result: Theorem 1.If P max,1 < P max,2 , then the positive equilibrium point (Q * , P * , E * ) is unstable with a two-dimensional stable manifold in parameter space.
If P max,1 = P max,2 , then the equilibrium point has the eigenvalue 0 and two negative eigenvalues.Thus there exists a one-dimensional center manifold which in this case consists of all equilibria, namely the equilibrium manifold M defined earlier.
For P max,1 > P max,2 , the coefficients a 0 , a 1 and a 2 are all positive.Furthermore This expression is positive if we make the following assumption (B): Note from Equations ( 2) and (3) that δ P E max,2 represents the maximal size of the immune effect E on P while δ Q E max,1 represents the maximal size of the immune effect E on Q.This effect is assumed to be stronger on the proliferating class of cells than on the quiescent class of cells.Thus assumption ( 21) is a natural one to make.This assumption is invariant under the actions of the drugs: so that these terms are multiplied by the same coefficients.Hence we also have the following result: Theorem 2. If P max,1 > P max,2 and δ P E max,2 ≥ δ Q E max,1 , then the positive equilibrium point (Q * , P * , E * ) is locally asymptotically stable.
The limiting case P max,1 = P max,2 represents a degenerate scenario.In many cases no positive equilibrium exists.For example, if Also, although the positive equilibrium point in Theorem 2 is stable, the value can be very high.In fact, P * diverges to +∞ as these parameter relations are reached (c.f. ( 18)): .
For the equilibrium values to be relatively small ('chronic'), we see that P max,1 must be significantly larger than P max,2 .In terms of the parameter values with drug actions, this can be achieved using the drugs u 2 and u 3 which increase Pmax,1 and decrease Pmax,2 , c.f., So drug administration shifts the balance towards P max,1 and this creates an asymptotically stable positive equilibrium point (Q * , P * , E * ), hopefully with low values for P * and Q * .
Figure 4 shows how the positive equilibrium values change as (only) one of the controls is varied.Note that the equilibrium values for Q and E do not change if only the control u 1 is varied.Also for changes in the controls u 2 and u 3 these equilibrium values change little and in the graphs the corresponding curves are almost constant.However, in these cases the equilibrium values for Q and P are well-controlled by the therapies.Contrary to the case when Q = 0, the u 2 and u 3 controls have strong effects by cutting down the influx of cells from the Q into the P compartment.All equilibria shown in these graphs satisfy the conditions of Theorem 2 and are locally asymptotically stable.  1 and 2.

QFigure 1 .
Figure 1.Diagram of the dynamical system.The green circular areas represent the "populations" included in the model.Solid arrows extending from or to the populations represent changes in numbers, with inward-pointing arrows representing increases and outward-pointing arrows decreases.Dashed arrows indicate indirect effects on those increases or decreases.Bars represent inhibition of a production or an indirect effect, due to the represented treatment; arrows represent amplification of a rate or an indirect effect.The effects of the general BCR-ABL1 inhibitor u 1 are shown using orange dashed bars and arrows, the effects of the BCR-ABL1 inhibitor u 2 which also has immune effects are shown using wide red solid bars and arrows and the effects of the immunomodulatory compound u 3 are shown using blue solid bars and arrows.

Figure 3 .
Figure 3. Phase portraits of the reduced dynamics for Q = 0 and P max,1 > P max,2 for the uncontrolled system (top, left) and for constant controls u 1 ≡ 1 (top, right), u 2 ≡ 1 (bottom, left) and u 3 ≡ 1 (bottom, right).The numerical values for these phase portraits are given in Tables1 and 2.

Figure 4 .
Figure 4.The values of the positive equilibrium point (Q * , P * , E * ) as the values for a single control are varied from 0 to 1.The parameter values used in the computations are given in Tables1 and 2.

Table 1 .
States and parameters for the dynamical system.

Table 2 .
Controls and pharmacodynamic parameters.
+ P max,2 ) ; then E ∆ ≤ E ss and for E < E ∆ we have that dE dt > 0. Furthermore, for E = E ss , Equation (8) reduces to dE dt = −δ E 1 + E < 0 and thus the values of E cannot reach the value E ss if they start below E ss .

Proposition 2 .
If P max,2 > P max,1 , then multiple equilibria (P * , E * ) inside P 0 can exist.At points (P saddle-node bifurcations occur which a stable and an unstable equilibrium point merge.(Pmax,2 − P max,1 ) δ E E * (1 + P * ) (1 + P * + P max,1 P * , E * ) where δ P r P (P max,2 − P max,1 ) P * (1 + P * + P max,1 P * ) 2 * (1 + E * ) 2 = 1 (14) * ) then a positive equilibrium point (Q * , P * , E * ) exists, but this relation is non-generic and generally will not be satisfied for a given set of parameters.3.5.Existence and Stability of a Positive Equilibrium Point (Q * , P * , E * ) for P max,1 = P max,2 If P max,1 < P max,2 , this solution is positive if and only if If one of these inequalities is violated, no positive equilibrium solution P * = P * (E * ) exists and the overall dynamics are determined either by exponential growth or decay of Q.If P * = P * (E * ) exists and is positive, then Equation (7) defines Q * ask P Q * = δ P 1 + E max,2 * P * .