Oncolytic virotherapy is a cancer treatment that uses viruses to selectively infect and destroy cancer cells while leaving healthy cells unharmed. These viruses, known as oncolytic viruses, replicate within tumor cells, causing them to lyse (burst) and release new viral particles that can infect surrounding cancer cells. This process also releases tumor antigens, which can trigger an immune response against the cancer. The dynamics of Oncolytic virotherapy are highly nonlinear. 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 oncolytic dynamic 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 a Hopf bifurcation point in one of the models and branch points in all the three models. The Hopf bifurcation point was eliminated using an activation factor that involves the tanh function.  The branch points (which cause multiple steady-state solutions from a singular point) are very beneficial because they enable the Multiobjective nonlinear model predictive control calculations to converge to the Utopia point (the best possible solution) in the models. It is proved (with computational validation) that the branch points were caused because of the existence of two distinct separable functions in one of the equations in each dynamic model. A theorem was developed to demonstrate this fact for any dynamic model.

Mullen et al (2002) [1] describe the mechanism of viral oncolysis, while Chiocca (2002) [2] present several oncolytic viruses, and Aghi et al (2005) [3] discuss the various oncolytic viral therapies. Novozhilov et al (2006) [4] mathematically modeled tumor therapy with oncolytic viruses. Kelly, and Russell (2007) [5] present the history of oncolytic viruses from genesis to genetic engineering. Alonso et al (2007) [6] showed that the combination of the oncolytic adenovirus ICOVIR-5 with chemotherapy provides an enhanced anti-glioma effect. Zurakowski, and Wodarz (2007) [7], discuss various model-driven approaches for in vitro combination therapy using ONYX-015 replicating oncolytic adenovirus. Wong et al (2010) [8] present strategies for overcoming the obstacles in using oncolytic viruses for cancer therapy. Ottolino et al (2010) [9] provide intelligent designs combining therapy with oncolytic viruses. Komarova and Wodarz (2010) [10] present ODE models for oncolytic virus dynamics. Eager and J. Nemunaitis (2011) [11] researched the various clinical development directions in oncolytic viral therapy. Agarwal and Bhadauria (2011[12]) modelled and analyzed tumor therapy with oncolytic virus. Bagheri et al (2011) [13] provide a dynamical systems model for combinatorial cancer therapy that enhances oncolytic adenovirus efficacy by MEK-inhibition. Tian (2011) [14] shows the replicability of oncolytic virus, defining conditions in tumor virotherapy. Donnelly et al (2012) [15] present the recent clinical experiences with oncolytic viruses. Russell et al (2012) [16] provide a discussion of oncolytic virotherapy. Zhou et al (2013) [17] present the clinical research progress for oncolytic adenovirus targeting cancer therapy. Patel and Kratzke (2013) [18] discuss the first wave of translational clinical trials in oncolytic virus therapy for cancer. Wang et al (2013) [19] show that the lytic cycle is a defining process in oncolytic virotherapy. Kim et al (2015) [20] discuss the quantitative impact of immunomodulation versus oncolysis with a cytokine-expressing virus therapeutic. Si and Zhang (2015) [21] discuss the control exponential growth of tumor cells with the slow spread of oncolytic virus. Chen and Su (2016) [22] provide an improved model of tumor therapy with oncolytic virus. Simpson et al (2016) [23] discussed the recent advances in cancer immunotherapy via combining oncolytic virotherapy with chemotherapy. Su et al (2016) [24] developed an optimal control model of tumor treatment with oncolytic virus and MEK Inhibitor.  Malinzi et al (2018) [25] developed a mathematical and optimal control analysis on the enhancement of chemotherapy using oncolytic virotherapy. Adi-Kusumo et al (2020) [26] showed the existence of a Hopf Bifurcation on a cancer therapy model by oncolytic virus involving the malignancy effect and therapeutic efficacy. This work aims to perform bifurcation analysis and multiobjective nonlinear control (MNLMPC) studies in three oncolytic virus models, which are discussed in Adi-Kusumo et al (2020) [26] (model 1) Su et al (2016) [24] (model 2), and Malinzi et al (2018) [25] (model 3). The paper is organized as follows. First, the model equations are presented, followed by a discussion of the numerical techniques involving bifurcation analysis and multiobjective nonlinear model predictive control (MNLMPC). The results are then presented, followed by the discussion and conclusions.

Model Equations 

Model 1 (Adi-Kusumo et al (2020) [26]) 

The ODE set representing the first model is 

xval and yval represent the uninfected cancer cells and the infected ones by oncolytic viruses. 

The base parameters are 

Model 2 Su et al (2016) [24] (model 2)

The model equations are

xval, yval, and vval stand for the population of uninfected cells, infected tumor cells, and oncolytic viruses, and zval represents the average expression level of CAR on the surface of the cells.

The base parameter values are

Model 3Malinzi et al (2018) [25]

Here, (uval, iva, vval,cval) represent the  un-infected tumour density, virus-infected tumour cell density, free virus particles, and, drug concentration respectively. u1 and u2 are the control parameters. 

The base parameter values 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[27]; Dhooge Govearts, Kuznetsov, Mestrom and Riet, 2004[28]). This program detects Limit points (LP), branch points (BP), and Hopf bifurcation points(H) for an ODE system.

  Let the bifurcation parameter be.  Since the gradient is orthogonal to the tangent vector, 

The tangent plane at any point      must satisfy 

Where A is 

where is the Jacobian matrix. For both limit and branch points, the matrix must be singular. The n+1 th component of the tangent vector = 0 for a limit point (LP)and for a branch point (BP) the matrix must be singular. At a Hopf bifurcation point, 

@ 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[29]; 2009 [30]) and Govaerts [2000] [31]. 

Hopf bifurcations cause unwanted oscillatory behavior and limit cycles. The tanh activation function (where a control value u is replaced by) is commonly used in neural nets (Dubey et al 2022[32]; Kamalov et al, 2021[33] and Szandała, 2020[34]) and optimal control problems (Sridhar 2023[35]) to eliminate spikes in the optimal control profile. Hopf bifurcation points cause oscillatory behavior. Oscillations are similar to spikes, and the results in Sridhar (2024) [36] demonstrate that the tanh factor also eliminates the Hopf bifurcation by preventing the occurrence of oscillations. Sridhar (2024) [36] explained with several examples how the activation factor involving the tanh function successfully eliminates the limit cycle causing Hopf bifurcation points. This was because the tanh function increases the time period of the oscillatory behavior, which occurs in the form of a limit cycle caused by Hopf bifurcations.

Multiobjective Nonlinear Model Predictive Control (MNLMPC) 

Flores Tlacuahuaz et al (2012) [37] developed a multiobjective nonlinear model predictive control (MNLMPC) method that is rigorous and does not involve weighting functions or additional constraints. This procedure is used for performing the MNLMPC calculations Here (j=1, 2.n) represents  the variables that need to be minimized/maximized simultaneously for a problem involving a set of ODE.

being the final time value, and n the total number of objective variables and. u the control parameter. This MNLMPC procedure first solves the single objective optimal control problem independently optimizing each of the variables individually. The minimization/maximization of will lead to the values . Then the optimization problem that will be solved is 

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 ( for all j) is obtained. 

Pyomo (Hart et al, 2017) [38] 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) [39]and confirmed as a global solution with BARON (Tawarmalani, M. and N. V. Sahinidis 2005) [40]. 

The steps of the algorithm are as follows

1. Optimize and obtain at various time intervals ti. The subscript i is the index for each time step. 
2. Minimize and get the control values for various times.
3. Implement the first obtained control values 
4. Repeat steps 1 to 3 until there is an insignificant difference between the implemented and the first obtained value of the control variables or if the Utopia point is achieved. The Utopia point is when for all j. 

 Sridhar (2024) [41] proved that the MNLMPC calculations to converge to the Utopia solution when the bifurcation analysis revealed the presence of limit and branch points. This was done by imposing the singularity condition on the co-state equation (Upreti, 2013) [42]. If the minimization of lead to the value and the minimization of lead to the value The MNLPMC calculations will minimize the function . The multiobjective optimal control problem is

Differentiating the objective function results in 

The Utopia point requires that both and are zero.  Hence  

the optimal control co-state equation (Upreti; 2013) [42] is 

is the Lagrangian multiplier. is the final time.  The first term in this equation is 0 and hence 

At a limit or a branch point, for the set of ODE is singular. Hence there are two different vectors-values for where . In between there is a vector where . This, coupled with the boundary condition will lead to This makes the problem an unconstrained optimization problem, and the only solution is the Utopia solution.  

Model 1

When pval was used as the bifurcation parameter a Hopf bifurcation point was located at (xval, yval, pval) values of (3.086053, 55.448945, -0.898380). This is shown in Fig. 1a. The limit cycle caused by this Hopf bifurcation point is shown in Fig. 1b.  When pval was changed to pval(tanh(pval)) the Hopf bifurcation point disappears, validating the analysis in Sridhar (2024) [36] where it is demonstrated with several examples that the activation factor involving the tanh function successfully eliminates the limit cycle causing Hopf bifurcation points.  This is shown in Fig. 1 c.  When r2 was used as the bifurcation parameter a branch point was located at 

(xval, yval, r2) values of (0.000000 50.049900 4.003996). This is shown in Figure. 1d. 

For the MNLMPC calculations, were    minimized individually and led to values of 200 and 0. pval was the control parameter. The multiobjective optimal control problem will involve the minimization of subject to the equations governing Model 1. This led to a value of zero (the Utopia solution).   The MNLMPC value of pval was 0.54427. Figs 1e, 1f and 1g. show the various MNLMPC profiles.

Model 2

In model 2, when bifurcation analysis was performed, with u as a parameter, a branch point was located at (xval, yval, vval,zval,u) values of ( 9.e+08, 0,0, 1.043561 0.582576 ). This is seen in Fig. 2a. 

For the MNLMPC calculations, xval(0) was set as 400 and    were minimized individually, yielding values of 400, 0 and 0. The multiobjective optimal control problem will involve the minimization of subject to the equations governing Model 2. This led to a value of zero (the Utopia solution).  The MNLMPC value of u was 0.932. Figs 2b,2c,2d, and 2e show the various MNLMPC profiles for Model 2.

Model 3

In model 3, when bifurcation analysis was performed, with  as a parameter, a branch point was located at values of. (0, 50, 0, 0.119904, 0.000065). Additionally, a limit point was found at values of (0.009703, 50.002421, 0, 0.119904, 0.000065). This is shown in Figure. 3a. 

For the MNLMPC calculations, uval(0) was set as 1000, ival(0) was set as 75.  were minimized individually, yielding values of 6.1539610796052817e+04, and 75. The multiobjective optimal control problem will involve the minimization of subject to the equations governing Model 3. This led to a value of zero (the Utopia solution).  The MNLMPC value of u1 and u2 (the two control parameters) were 7.863627108510299e-07 and 16.542884153050405. Figs 3b-3g show the various MNLMPC profiles for Model 3.

                                                                           Figure 1a: Bifurcation Analysis in Model 1 indicating Hopf bifurcation

                                                                Figure 1b: Bifurcation Analysis in Model 1 indicating limit cycle caused by Hopf bifurcation

                                                   Figure 1c: Bifurcation Analysis in Model 1 (Hopf bifurcation eliminated by tanh activation function)

                                                                     Figure 1d: Bifurcation Analysis in Model 1 indicating Branch Point

                                                                                                   Figure 1e: MNLMC (Model 1) yval vs t

                                                                                                   Figure 1f: MNLMC (Model 1) xval vs t

                                                                                                         Figure 1g: MNLMC (Model 1) pval vs t

                                                                        Figure 2a: Bifurcation Analysis in Model 2 indicating branch point

                                                                                                             Figure. 2b: MNLMC (Model 2) xval vs t

                                                                                                       Figure 2c: MNLMC (Model 2) yval vs t

                                                                                                        Figure 2d: MNLMC (Model 2) vval vs t

                                                                                                      Figure 2e: MNLMC (Model 2) u vs t

                                                                 Figure 3a: Bifurcation Analysis in Model 3 indicating Branch Point and limit point

                                                                                                       Figure 3b: MNLMPC in Model 3 vval vs t

                                                                                                         Figure 3c: MNLMPC in Model 3 ival vs t

                                                                                                 Figure 3d: MNLMPC in Model 3 cval vs t

                                                                                                           Figure 3e: MNLMPC in Model 3 uval vs t

                                                                                                    Figure 3f: MNLMPC in Model 3 u1 vs t

                                                                                                    Figure 3g: MNLMPC in Model 3 u2 vs t

Model 1 exhibits a Hopf bifurcation point, which is eliminated using an activation factor involving a tan function, validating the analysis in Sridhar (2024) [36]. All three models exhibit branch points. These branch points enable the MNLMC calculations to yield the Utopia solution, validating the analysis in Sridhar (2024) [41]. The following paragraphs explain the cause of the occurrence of the branch points. 

Theorem

If one of the functions in a dynamic system is separable into two distinct functions, a branch point singularity will occur in the system. 

Proof

Consider a system of equations

 . Defining the matrix A as 

is the bifurcation parameter. The matrix A can be written in a compact form as 

The tangent at any point x; must satisfy 

The matrix must be singular at both limit and branch points. The n+1 th component of the tangent vector = 0 at a limit point (LP) and for a branch point (BP) the matrix must be singular. 

Let any of the functions fi are separable into 2 functions as 

At steady-state and this will imply that either or or both and  must be 0.  This implies that two branches and  will meet at a point where both and are 0. 

At this point, the matrix B will be singular as a row in this matrix would be  

The singularity in B implies that there exists a branch point.

In model 1, a branch point was located at 

(xval, yval) values of (0.000000 50.049900). Here, the two distinct functions can be obtained from the first equation in Model 1 

The 2 distinct function branches are

and 

The values satisfy both the above equations confirming the correctness of the theorem.

 In model 2, a branch point was located at (xval, yval, vval,zval,u) values of ( 9.e+08; 0;0; 1.043561 0.582576 ). The two distinct functions can be seen from the second equation of model 2, 

Here, the two distinct function branches are 

and 

Substituting  

Satisfies both the above equations, confirming the correctness of the theorem.

In model 3, a branch point was located at values of. (0, 50, 0, 0.119904, 0.000065) and the two distinct branches can be obtained from the first equation in model 3 which is 

These 2 branches are 

and

Substituting

satisfies both the above equations confirming the correctness of the theorem.

Bifurcation analysis and multiobjective nonlinear control (MNLMPC) studies in three oncolytic virus models. The bifurcation analysis revealed the existence of a Hopf bifurcation point in one of the models and branch points in all the three models. The Hopf bifurcation point was eliminated using an activation factor that involves the tanh function.  The branch points (which cause multiple steady-state solutions from a singular point) are very beneficial because they enable the Multiobjective nonlinear model predictive control calculations to converge to the Utopia point (the best possible solution) in the models. It is proved (with computational validation) that the branch points were caused because of the existence of two distinct separable functions in one of the equations in each dynamic model. A theorem was developed to demonstrate this fact for any dynamic model.  A combination of bifurcation analysis and Multiobjective Nonlinear Model Predictive Control (MNLMPC) for dynamic models involving oncolytic viral therapy is the main contribution of this paper. 

All data used is presented in the paper

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

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