Millions of people are affected by leukemia. It is important to understand the progression dynamics of this disease to be able to minimize the damage that is caused by it. This article provides a mathematical framework to develop strategies to control leukemia. Bifurcation analysis is a powerful mathematical tool used to deal with the nonlinear dynamics of any process. Several factors must be considered, and multiple objectives must be met simultaneously.  Bifurcation analysis and multiobjective nonlinear model predictive control (MNLMPC) calculations are performed on three leukemia models. The MATLAB program MATCONT was used to perform the bifurcation analysis. The MNLMPC calculations were performed using the optimization language PYOMO in conjunction with the state-of-the-art global optimization solvers IPOPT and  BARON. The bifurcation analysis revealed the existence of limit points and branch in the models. The limit and branch points were beneficial because they enabled the multiobjective nonlinear model predictive control calculations to converge to the Utopia point in both problems, which is the most beneficial solution. A combination of bifurcation analysis and multiobjective nonlinear model predictive control for leukemia models is the main contribution of this paper. 

Deininger et al (2000)[1] investigated the molecular biology of chronic myeloid leukemia.  Topaly et al [2] researched the synergistic activity of the new ANL-specific tyrosine kinase inhibitor STI571 and the effect of chemotherapeutic drugs on BCR-ABL-positive chronic myelogenous leukemia cells. Deininger and co-workers(2003)[3,4] studied the effect of imatinib on patients with chronic myeloid leukemia. Goldman and Melo(2003) [5] described the treatments for chronic myeloid leukemia patients.  Mackey et al (2004)[6] studied the periodic behavior of chronic myelogenous leukemia. Baccarani et al [7] studied the effect of  Imatinib and pegylated human recombinant interferon-alpha2 b on early chronic-phase chronic myeloid leukemia. Moore and  Li(2004)[8] developed a mathematical model of chronic myelogenous leukemia. Adimy et al (2005)[9] performed a mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia.  Michor et al (2005)[10], studied the dynamics of chronic myeloid leukemia. Nanda et al (2007)[11] , performed single objective optimal control of a mathematical model of chronic myelogenous leukemia.   Liso et al and Ommen et al (2008)[12,13] and Cucuianu et al (2010)[14] performed theoretical studies of mathematical models involving myelogenous leukemia. Dohner et al (2010)[15] provided recommendations for the diagnosis and management of acute myeloid leukemia in adults.  Komarova(2011)[16], Stiehl(2012)[17], Maclean et al (2013,2014) [18,19], Agarwal (2015)[20], and Clapp(2015)[21] performed mathematical investigations of problems involving leukaemia.    Crowell et al (2016)[22] studied feedback mechanisms control coexistence in a stem cell model of acute myeloid leukaemia.  Austin et al (2016)[23] , described harnessing the immune system in acute myeloid leukaemia. Zeidan et al (2016)[24] described the economic burden associated with acute myeloid leukemia treatment.Masarova et al (2017) [25] performed additional investigations about harnessing the immune system against leukemia. Lichtenegger et al (2017) [26] presented more developments in immunotherapy of acute myeloid leukemia.  Krupar et al (2018) [27], provided an   analysis of anti-leukemia immune response and  immune evasion in acute myeloid leukemia. Sharp et al (2019)[28] and Khatun et al (2020) [29] performed single-objewctive optimal control studies of ,myeloid leukaemia treatment , Journal of Theoretical Biology, Volume 470, 2019, Pages 30-42, ISSN 0022-5193. In this work, bifurcation analysis is performed in conjunction with multiobjective nonlinear model predictive control (MNLMPC) for three leukemia models that are described in Nanda et al (2007)[11] , Khatun et al (2020), and Sharp et al (2019) (Model 1, Model 2, and Model 3).   This paper is organized as follows. First, the leukemia  models   are presented.  The numerical procedures (bifurcation analysis and multiobjective nonlinear model predictive control (MNLMPC) are then described. This is followed by the results and discussion, and conclusions. 

Model 1

The model equations are 

Bifurcation analysis 

The MATLAB software MATCONT is used to perform the bifurcation calculations. Bifurcation analysis deals with multiple steady-states and limit cycles.  Multiple steady states occur because of the existence of branch and limit points.  Hopf bifurcation points cause limit cycles.  A commonly used MATLAB program that locates limit points, branch points, and Hopf bifurcation points is MATCONT (Dhooge Govearts, and Kuznetsov, 2003[30]; Dhooge Govearts, Kuznetsov, Mestrom and   Riet, 2004[31]).  This program detects Limit points (LP), branch points (BP), and Hopf bifurcation points(H) for an ODE system

where @ Indicates the bialternate product while  is the n-square identity matrix. Hopf bifurcations cause limit cycles and should be eliminated because limit cycles make optimization and control tasks very difficult.  More details can be found in Kuznetsov (1998[32]; 2009[33]) and Govaerts [2000] [34]

Multiobjective Nonlinear Model Predictive Control (MNLMPC) 

This will provide the values of u at various times. The first obtained control value of u is implemented and the rest are discarded. This procedure is repeated until the implemented and the first obtained control values are the same or if the Utopia point where is obtained.

Pyomo (Hart et al., 2017) [36] is used for these calculations.  Here, the differential equations are converted to a Nonlinear Program (NLP) using the orthogonal collocation method   The NLP is solved using IPOPT (Wächter And Biegler, 2006) [37]and confirmed as a global solution with BARON (Tawarmalani, M. and N. V. Sahinidis 2005) [38]. 

The steps of the algorithm are as follows

For the bifurcation analysis of model 1, u1 is the bifurcation parameter and the other  parameter values are  sn=0.29; dn=0.35; de=0.40; dc=0.012;

This is shown in Fig. 1a. 

When u2 is the bifurcation parameter and the other parameter values are 

Figure 1a: (Bifurcation analysis model 1 u1 is bifurcation param

eter).

Figure 1b: (Bifurcation analysis model 1 u2 is bifurcation parameter).

For the bifurcation analysis of model 2, when u1 is the bifurcation parameter and the other parameter values are 

A branch point occurred at [sval, ival, wval, u1] values of (6.200000 0.000000 0.666667 0.003226). This is shown in Figure 2a.

Figure 2a: (Bifurcation analysis model 2 u1 is bifurcation parameter).

When u2 is the bifurcation parameter and the other parameter values are 

A branch point occurred at [sval, ival, wval, u2] values of (8.200000 0.0 0.0 0.0010). This is shown in Figure. 2b.

Figure 2b: (Bifurcation analysis model 2 u2 is bifurcation parameter).

For the bifurcation analysis of model 3, u is the bifurcation parameter and the other parameter values are 

This is shown in Figure.3.

Figure 3: (Bifurcation analysis model 3 u is bifurcation parameter).

Figure 4a:(MNLMPC for model 1 tn vs t)

Figure 4b: (MNLMPC for model 1 tn vs t)

Figure 4c: (MNLMPC for model 1 c vs t)

The obtained control profile of u1 and u2 exhibited noise (Figure. 4d and Figure. 4e). This issue was addressed using the Savitzky-Golay Filter. The smoothed version of this profile is shown in Figure 4f and 4g. The MNLMPC control values obtained for u1 and u2 are 0.0624 and 0.00443. The MNLMPC calculations converged to the Utopia solution, validating the analysis by Sridhar (2024), which demonstrated that the presence of a limit point/branch point enables the MNLMPC calculations to reach the optimal (Utopia) solution.

Figure 4d: (MNLMPC for model 1 u1 vs t)

Figure 4e:(MNLMPC for model 1 u2 vs t)

Figure 4f:(MNLMPC for model 1 u1 (Savitzky Golay) vs t)

Figure 4g: MNLMPC for model 1 u2  (Savitzky Golay)  vs t)

subject to the equations governing the model. This led to a value of zero (the Utopia solution. The various concentration profiles for this MNLMPC calculation are shown in Figure 5a-5c. 

Figure 5a: (MNLMPC model 2, sval vs t)

Figure 5b: (MNLMPC model 2, ival vs t)

Figure 5c: (MNLMPC model 2, wval s t)

The obtained control profile of u1 and u2 exhibited noise (Figure. 5d and Figure. 5e). This issue was addressed using the Savitzky-Golay Filter. The smoothed version of this profile is shown in Figures 4f and 4g. The MNLMPC control values obtained for u1 and u2 are 0.002779 and 0.02749. The MNLMPC calculations converged to the Utopia solution, validating the analysis by Sridhar (2024), which demonstrated that the presence of a limit point/branch point enables the MNLMPC calculations to reach the optimal (Utopia) solution.

Figure 5d: (MNLMPC model 2, sval vs t)

Figure 5e:(MNLMPC model 2, u2 vs t)

was minimized subject to the equations governing the model. This led to a value of zero (the Utopia solution. The MNLMPC control values obtained for u was 0.6345 The various concentration profiles for this MNLMPC calculation are shown in Figs. 6a-6f. The MNLMPC calculations converged to the Utopia solution, validating the analysis by Sridhar (2024) [39], which demonstrated that the presence of a limit point/branch point enables the MNLMPC calculations to reach the optimal (Utopia) solution. 

Figure 5f: (MNLMPC model 2, u1 (Savitzky Golay) vs t)

Figure 5g: (MNLMPC model 2, u2 (Savitzky Golay) vs t)

Figure 6a: (MNLMPC model 3, sval vs t)

Figure 6b:(MNLMPC model 3, aval vs t)

Figure 6c:(MNLMPC model 3, tval vs t)

Figure 6d: (MNLMPC model 3, lval vs t)

Figure 6e: (MNLMPC model 3, dval vs t)

Figure 6f: (MNLMPC model 3, u vs t)

Model 1 and Model 3 displayed limit points, while Model 2 demonstrated a branch point. In all three instances, the MNLMPC calculations converged to the Utopia solution, validating the analysis of Sridhar (2024), which showed that the existence of a limit point or a branch point allows the MNLMPC calculations to achieve the best possible (Utopia) solution. 

Bifurcation analysis and Multiobjective nonlinear model predictive control calculations were performed on three leukemia models.   The bifurcation analysis revealed the existence of limit and branch points.  The limit and branch points (which cause multiple steady-state solutions from a singular point) are very beneficial as they enable the Multiobjective nonlinear model predictive control calculations to converge to the Utopia point (the best possible solution) in both models.  A combination of bifurcation analysis and multiobjective nonlinear model predictive control for leukemia disease models is the main contribution of this paper. 

Data Availability Statement

All data used is presented in the paper

Conflict of interest 

The author, Dr. Lakshmi N Sridhar has no conflict of interest.

Acknowledgement 

Dr. Sridhar thanks Dr. Carlos Ramirez and Dr. Suleiman for encouraging him to write single-author papers.