Bifurcation Analysis and Optimal Control of the Tumor Macrophage Interactions Volume 53- Issue 5

Lakshmi N Sridhar*

• Chemical Engineering Department, University of Puerto Rico, USA

Received: November 15, 2023;   Published: November 24, 2023

*Corresponding author: Lakshmi N Sridhar, Chemical Engineering Department, University of Puerto Rico, Mayaguez, PR 00681, USA

DOI: 10.26717/BJSTR.2023.53.008470

Abstract PDF

Bifurcation analysis and optimal control was performed on a model describing the interaction between tumors and macrophages. The MATLAB program MATCONT was used to perform the bifurcation analysis and the optimization language promo is used in conjunction with the state-of-the-art global optimization solvers IPOPT and BARON for the optimal control calculations. Both single and multi-objective nonlinear model predictive control (that involves Multi objective optimal control) calculations were performed. The bifurcation analysis revealed oscillation causing Hopf bifurcations while the optimal control showed the existence of spikes in the control profiles. Both were eliminated using an activation factor involving the hyperbolic tangent function.

Keywords: Bifurcation; Optimization; Control; Macrophage; Tumor

Genetically engineered macrophages that kill tumor cells cause inflammation, leading to DNA and tissue damage. This is remedied by alternatively activated macrophages that can reduce inflammation and repair tissue damage. The interaction between the tumor cells and the two macrophages demonstrates oscillatory behavior that from a mathematical standpoint arises from the existence of a Hopf bifurcation. Oscillations are similar to the existence of spikes in a control profile. Traditionally spikes in optimal control problems are eliminated using a Tanh activation factor. This raises the question of whether using the tanh activation factor will eliminate the unwanted oscillatory behavior that exists because of a Hopf bifurcation. In this work, both bifurcation analysis and optimal control are performed on a scaled dynamic model involving tumor cells and macrophages and it is shown that just as in optimal control, the tanh factor is effective in the elimination of spikes, it also eliminates the Hopf bifurcation that causes oscillatory behavior.

Mantovani and co-workers [1-3] studied tumor immunity, macrophage polarization, and the effect of Macrophage immunity on cancer. Wilke and Hahnfeldtt [4] developed a model studying the dichotomy of the immune response to cancer: Rakoff-Nahoum [5] studied the interaction between cancer and inflammation Prehn [6] studied investigated the effect of the immune reaction on the as a stimulator of tumor growth. Sica and co-workers [7] investigated Macrophage polarization in tumor progression. The two macrophages’ phenotypes are divided into M1 and M2 macrophages. M1 cells have anti-tumor properties that can produce pro-inflammatory cytokines to eradicate pathogens but can cause tissue and DNA damage. M2 cells can reduce the pro-inflammatory response and stimulate tissue repair. Allavena and Mantovani [8] suggested a re-polarization of macrophages towards the classically activated M1 cells, as an effective treatment approach to ensure tumor elimination. A considerable amount of computational work has been done dealing with tumor macrophage interactions. As early as 1985, De Boer, et al. [9] presented a model of the interactions between macrophages and T lymphocytes that generate an anti-tumor immune response. Owen and Sherratt [10] developed a five-dimensional differential equations model to study the roles of macrophages presence, influx, and ability to selectively kill tumor cells in avascular tumors. They verified that the proportion of macrophages in tumors increases with the chemotactic activity of the mutant cell line. Byrne, et al. [11] proposed a mathematical model that describes the interactions between normal cells, tumor cells and infiltrating macrophages, to evaluate the ability of the engineered macrophages to eliminate the tumor changes as model parameters vary Owen, et al. [12] investigated the role of chemotaxis, chemokine production and the efficacy of macrophages as vehicles for drug delivery to hypoxic tumor sites.

Webb, et al. [13] extended the model of Owen [12] to show that limited-diffusivity or non-cell- cycle dependent drugs help macrophages effectively target hypoxic tumor cells. Eftimie [14-16] studied the complex dynamics in the interactions between the macrophages and tumor cells. Dong, et al. [17] and Ghosh [18] and Liu, et al. [19] show that the models describing the dynamics of macrophage and tumor interactions exhibit periodic oscillations which indicate the presence of Hopf bifurcations. Shu, et al. [20] show the existence of Hopf bifurcations in the model describing the interaction between tumor cells, M1 and M2 macrophages. Sridhar [21] showed that the use of the hyperbolic tangent activation factor was effective in removing spikes in a control profile. Repeated spikes are similar to oscillatory behavior and the aim of this paper is to use a similar hyperbolic tangent activation factor on the tumor macrophage interaction model and demonstrate that the oscillation causing Hopf bifurcation can be eliminated. This paper is organized as follows. First, the tumor macrophage interaction model Shu, et al. [20] is described. This is followed by a discussion of the bifurcation analysis and the multi-objective nonlinear model predictive control strategy that involves optimal control. The results are then presented and discussed followed by conclusions.

Tumor Macrophage Interaction Model

The model and its scaled version Shu, et al. [20] is presented in this section. The model includes three variables: the population of tumor cells (T), the population of M1 macrophages (M1) and the population of M2 macrophages (M2). The interactions among these three cells is described by the following differential equations:

Owing to the large variation of parameters convergence for bifurcation analysis and optimal control is difficult. Hence a scaled model

Bifurcation Analysis

There has been a lot of work in chemical engineering involving bifurcation analysis throughout the years, The existence of multiple steady-states and oscillatory behavior in chemical processes has led to a lot of computational and analytical work to explain the causes for these nonlinear phenomena. Multiple steady states are caused by the existence of branch and limit points while oscillatory behavior is caused by the existence of Hopf bifurcations points. In the case of branch points and limit points the Jacobian matrix of the set of steadystate equations is singular. However, at a branch point, there are 2 distinct tangents at the singular point while at a limit point, there is only one tangent at the singular point. Singularities in the Jacobian matrix is often indicative of an optimal solution, and this motivates the investigation of how the singular points in the Jacobian matrix, indicated by branch and limit points would affect the multi-objective dynamic optimization. One of the most used software to locate limit points, branch points, and Hopf bifurcation points is MATCONT Dhooge, et al. [22,23]. This software detects Limit points, branch points and Hopfbifurcation points. Consider an ODE system:

for both limit and branch point the matrix B must be singular. For a limit point (LP) the n +1^th component of the tangent vector 1 0 n v + = and for a branch point (BP) the matrix must be singular. For a Hopf bifurcation, the function should be zero. indicates the biternate product while n I is the n-square identity matrix. More details can be found in Kuznetsov [24,25] and Govaerts [26]. The PI has published several articles where bifurcation analysis was used in chemical engineering problems such as reactive distillation and fermentation Ruiz, et al. [27-30]. Hopf bifurcation points cause unwanted oscillatory behavior which makes control tasks difficult.

Multi-Objective Nonlinear Model Predictive Control (MNLMPC) Method

The multi-objective nonlinear model predictive control (MNLMPC) method was first proposed by Flores Tlacuahuaz, et al. [31] and used by Sridhar [32-36]. This method is rigorous, and it does not involve the use of weighting functions do not do it impose additional parameters or additional constraints on the problem unlike the weighted function or the epsilon correction method (Miettinen; 1999). For a problem that is posed as

The MNLMPC method first solves dynamic optimization problems independently minimizing/maximizing each x_i individually. The minimization/maximization of x_i will lead to the values * x_i Then the optimization problem that will be solved is:

This is the MOOC calculation. It will provide the control values for various times. The first obtained control value is implemented and the remaining discarded. This procedure is repeated until the implemented and the first obtained control value are the same. The optimization package in Python, Pyomo Hart, et al. [37] where the differential equations are automatically converted to a Nonlinear Program (NLP) using the orthogonal collocation method Biegler, [38] is commonly used for these calculations. The state-of-the-art solvers like IPOPT Wachter and Biegler [39] and BARON Tawaralmani and Sahinidis [40] are normally used in conjunction with PYOMO. To summarize the steps of the algorithm are as follows.

1. Minimize/maximize i x subject to the differential and algebraic equations that govern the process using Pyomo and Baron. This will lead to the value x_i^*
2. Minimizemin (multi-objective function) subject to the differential and algebraic equations that govern the process. This is the MOOC calculation and provides the control values for various times. If this calculation results in obtaining a value of zero for the multi-objective function, then the Utopia point is obtained, and the cal
culations are terminated. Otherwise, we proceed to step 3. 3. Implement the first obtained control values and discard the remaining.
4. Repeat steps 1 to 3 until there is an insignificant difference between the implemented and the first obtained value of the control variables.
5. This strategy was used in Sridhar [30-34].

First, the bifurcation analysis was performed on a scaled model describing the tumor macrophage interactions. As shown by Shu [20], the model exhibits Hopf bifurcations when ζ is the bifurcation parameter. However, when the activation factor tanh (ζ ) ε is introduced and the resulting bifurcation parameter is ζ tanh (ζ )/ ε the Hopf bifurcation point disappears. ε is taken as 10−3 .Such an activation factor is commonly used to control spikes or miniature waves in optimal control problems. Figures 1a & 1b show the two profiles with and without the activation factor tanh (ζ ) ε . It can be seen that the use of the activation factor eliminated the Hopf bifurcation point. Figures 2a & 2b show the optimal control profiles for the single objective optimal control is performed minimizing ζ versus time minimizing the scaled variable X subject to the equations 4 5 and 6. The resulting value of the objective function is 8e-05. However, the control profile exhibits a spike as seen in Figure 2a. When the activation factor is used the spikes disappear as seen in Figure 1b. Using this corrected value of ζ Oscillatory behavior is like the existence of spikes, and it is shown that the activation factor involving the tanh function is effective in eliminating both and this is the main message in this article. The Multi-Objective Nonlinear Model Predictive Control (MNLMPC) calculations were also performed. First, X was minimized while Y and Z were maximized individually [41-43]. The objective function values obtained are 0.0109, 8.95 and 2.156. The objective function for the MNLMPC calculations is subject to the scaled differential equations. The first obtained control value of value of x was implemented and the process is repeated until there is no difference between the first obtained control value and the implemented control value. In this case tis value is 0.0292. Figures 3-6 show the profiles of the various variables where MNLMPC was used.

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Bifurcation analysis and optimal control was performed on a scaled model involving interactions between macrophages and tumors. Oscillation causing Hopf bifurcation points were found. It is shown that the incorporation of the hyperbolic tangent function (Tanh) that is normally used to eliminate spikes in optimal control profiles is also effective in the elimination of the unwanted oscillation causing Hopf bifurcation points.

All data used is presented on paper.

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

1. Mantovani A, Romero P, Palucka AK, Marincola (2008) Tumour immunity: effector response to tumour and role of the Lancet 371(9614): 771-783.
2. Mantovani A, Sozzani S, Locati M, Allavena P, Sica A (2002) Macrophage polarization: Tumor-associated macrophages as a paradigm for polarized M2 mononuclear phagocytes. Trends Immunol 23(11): 549-555.
3. Mantovani A, Sica A (2010) Macrophages, innate immunity and cancer: balance, tolerance, and diversity. Curr Opin Immunol 22(2): 231-237.
4. Wilkie KP, Hahnfeldt P (2017) Modeling the dichotomy of the immune response to cancer: cytotoxic effects and tumor-promoting Bull Math Biol 79(6): 1426-1448.
5. Rakoff Nahoum S (2006) Why cancer and inflammation? Yale J Biol Med 79: 123-130.
6. Prehn RT (1972) The immune reaction as a stimulator of tumor Science 176(4031): 170-171.
7. Sica A, Larghi P, Mancino A, Rubino L, Porta C, et al. (2008) Macrophage polarization in tumour Semin Cancer Biol 18(5): 349-355.
8. Allavena P, Mantovani A (2012) Immunology in the clinic review series; Focus on cancer: Tumour-associated macrophages: Undisputed stars of the inflammatory tumour microenvironment. Clin Exp Immunol 167(2): 195-205.
9. De Boer RJ, Hogeweg P, Dullens HF, RAD Weger, Otter WD, et al. (1985) Macrophage T lymphocyte interactions in the anti-tumor immune response: A mathematical model. J Immunol 134(4): 2748-2758.
10. Owen MR, Sherratt JA (1998) Modelling the macrophage invasion of tumours: Effects on growth and IMA J Math Appl Med Biol 15(2): 165-185.
11. Byrne HM, Cox SM, Kelly CF (2004) Macrophage-tumour interactions: In vivo Discrete Contin. Dyn Syst Ser B 4(1): 81-98.
12. Owen MR, Byrne HM, Lewis CF (2004) Mathematical modelling of the use of macrophages as vehicles for drug delivery to hypoxic tumour J Theor Biol 226(4): 377-391.
13. Webb SD, Owen MR, Byrne HM, Murdoch C, Lewis CF (2007) Macrophage- based anti-cancer therapy: Modelling different modes of tumour targeting. Bull Math Biol 69(5): 1747-1776
14. Den Breems NY, Eftimie R (2016) The re-polarisation of M2 and M1 macrophages and its role on cancer outcomes. J Theor Biol 390: 23-39.
15. Eftimie R, Hamam H (2017) T cells macrophages paradox in melanoma immunotherapies. J Theor Biol 420: 82-104.
16. Eftimie R, Eftimie G (2018) Tumour-associated macrophages and oncolytic virotherapies: A mathematical investigation into a complex Lett Biomath 5(1): 70-99.
17. Dong Y, R Miyazaki, G Huang, Y Takeuchi (2015) Dynamics in a tumor immune system with time delays. Appl Math Comput 252: 99-113.
18. Khajanchi S, Perc M, Ghosh D (2018) The influence of time delay in a chaotic cancer model Chaos 28: 103101.
19. Liu D, Ruan S, Zhu D (2012) Stable periodic oscillations in a two-stage cancer model of tumor and immune system interactions. Math Biosci Eng 9(2): 347-368.
20. Shu Yaqin, Jicai Huang, Yueping Dong, Yasuhiro Takeuchi (2020) Mathematical modeling and bifurcation analysis of pro- and anti-tumor Applied Mathematical Modelling 88: 758-773.
21. Sridhar LN (2023) Multi Objective Nonlinear Model Predictive Control of Diabetes Models Considering the Effects of Insulin and Archives Clin Med Microbiol 2(2): 23-32.
22. Dhooge A, Govearts W, Kuznetsov AY (2003) MATCONT: A Matlab package for numerical bifurcation analysis of ODEs. ACM transactions on Mathematical software 29(2): 141-164.
23. Dhooge A, Govaerts W, Kuznetsov YA, Mestrom W, Riet AM, et al. (2004) CL_MATCONT; A continuation toolbox in Matlab.
24. Kuznetsov YA (1998) Elements of applied bifurcation Springer NY.
25. Kuznetsov YA (2009) Five lectures on numerical bifurcation Utrecht University NL.
26. Govaerts wJF (2000) Numerical Methods for bifurcations of dynamical SIAM.
27. Ruiz Gerardo, Lakshmi N Sridhar, Raghunathan Rengaswamy (2006) The Isothermal Isobaric Reactive Flash By published in I & EC Res 45: 6548-6554.
28. Ruiz Gerardo, Misael Diaz, Lakshmi N Sridhar (2008) Singularities in reactive separation processes. by I & EC Res 47: 2808-2816.
29. Sridhar LN (2011) Elimination of oscillations in fermentation AIChE J 57(9): 2397-2405.
30. Sridhar LN (2012) Using Magnetic Nanoparticles to Eliminate Oscillations in Saccharomyces cerevisiae Fermentation Processes. JSBS 2: 27-32.
31. Flores Tlacuahuac A (2012) Pilar Morales and martin riveral toledo; Multiobjective Nonlinear model predictive control of a class of chemical reactors. I & EC research 51: 5891-5899.
32. Sridhar LN (2019) Multi objective optimization and nonlinear model predictive control of the continuous fermentation process involving Saccharomyces Biofuels 13(2).
33. Sridhar LN (2020) Multi objective nonlinear model predictive control of pharmaceutical batch crystallizers, Drug Development, and Industrial Pharmacy 46: 2089-2097.
34. Sridhar LN (2021) Single and Multi-objective optimal control of epidemic models involving vaccination and treatment accepted for publication in journal of biostatistics and epidemiology. J Biostat Epidemiol 7(1): 25-35
35. Sridhar LN (2022) Single and multi-objective optimal control of the wastewater treatment process, Transactions of the Indian National Academy of Engineering 7: 1339-1346.
36. Sridhar LN (2023) Multi Objective Nonlinear Model Predictive Control of Diabetes Models Considering the Effects of Insulin and Exercise. Archives Clin Med Microbiol 2(3): 60-39.
37. Hart William E, Carl D Laird, Jean Paul Watson, David L Woodruff, Gabriel A Hackebeil, et al. (2017) Pyomo – Optimization Modeling in Python. Second Edition Springer, p. 67.
38. Biegler LT (2007) An overview of simultaneous strategies for dynamic Chem Eng Process Process Intensif 46: 1043-1053.
39. Wächter A, Biegler L (2006) On the implementation of an interior-point filter line- search algorithm for large-scale nonlinear Math Program 106: 25-57.
40. Tawarmalani M, NV Sahinidis (2005) A polyhedral branch-and-cut approach to global optimization. Mathematical Programming 103(2): 225-249.
41. Ghosh D, S Khajanchi, S Mangiarotti, F Denis, SK Dana, et al. (2017) How tumor growth can be influenced by delayed interactions between cancer cells and the microenvironment. BioSystems 158: 17-30.
42. Khajanchi S, Ghosh D (2015) The combined effects of optimal control in cancer remission. Appl Math Comput 271: 375-388.
43. Miettinen Kaisa (1999) Nonlinear Multiobjective Optimization; Kluwers international series.