Stability analysis of the singular points and Hopf bifurcations of a tumor growth control model

We carry out qualitative analysis of a fourth‐order tumor growth control model using ordinary differential equations. We show that the system has one positive equilibrium point, and its stability is independent of the feedback gain. Using a Lyapunov function method, we prove that there exist realistic parameter values for which the systems admit limit cycle oscillations due to a supercritical Hopf bifurcation. The time evolution of the state variables is also represented.


INTRODUCTION
Personalized therapies are gaining more interest in clinical practice due to their many advantages, like greater efficiency and fewer side effects.An efficient approach is to use mathematical models to describe specific physiological problems [1][2][3][4][5] and use control theory to generate therapies [6][7][8][9][10].Mathematical optimization and personalization of cancer chemotherapies may also decrease the chance of drug resistance development during the treatment [11].We consider the control problem of a tumor growth model and carry out qualitative analysis of the closed-loop system.
The literature of tumor models is rich; see, for example, the reviews by Altrock et al. [12] and Jarrett et al. [13].In our analysis, we use a tumor model describing the tumor growth with linear dynamics, using nonlinear pharmacodynamics and linear pharmacokinetics [14], discussed in Section 2. The model is formulated using formal reaction kinetics [15,16].We chose this model because it describes all the fundamental physiological phenomena in a compact model, and it has been successfully used for therapy optimization in mice experiments [17,18].
In our earlier works, we analyzed the closed-loop control of a simplified model [19], using state feedback [20].The qualitative analysis of the closed-loop model showed that there is one positive equilibrium.We also found inequalities in the design parameters to ensure the stability of the closed-loop system, that is, the efficacy of the therapy.However, there were no bifurcations in the closed-loop system.Later, we analyzed the closed-loop system based on a more complex model [14], with simple feedback based on the measured variable [21].We found one positive singular point and bifurcations for realistic parameter values.We continue the analysis further in Sections 3-5.
The tumor model and the control problem used for the analysis are described in Section 2. In Section 3, we find one singular point and realistic model parameter values, which result in pure imaginary eigenvalues of the Jacobian at the singular point.We also find that the eigenvalues of the Jacobian do not depend on the control gain.In Section 4, we analyze the stability of the equilibrium by the Lyapunov direct method [22].We examine the evolution of the states in Section 5 and draw the conclusions in Section 6.
In our analysis, we use the Wolfram Mathematica 13.0 computer algebra system for symbolic calculations, numerical verifications, and also for creating the figures.The Mathematica notebook is available online [23].

TUMOR GROWTH MODEL
We consider the control of tumor growth described by four ordinary differential equations [14].The four equations define the dynamics of the four states, that is, the x 1 time function of proliferating tumor volume, the x 2 time function of dead tumor volume, and x 3 and x 4 being the states of a two-compartment model.The state of the first compartment x 3 is the drug level in the central compartment (the blood), and the state of the second compartment is the drug level in the peripheral compartment x 4 .The model was created with an analogy to formal reaction kinetics [14][15][16]24], which improves communication and validation with medical experts.
The system is defined by the equations ( The measured quantity in practice is  = x 1 + x 2 , which is considered to be the output of the system.The practical input is drug injection/infusion into the blood (x 3 ), with injection rate u [mg•kg −1 •day −1 ].The tumor volumes have the dimension mm 3 , and the drug levels have the dimension mg•kg −1 .Table 1 shows a possible value of model parameters, whose values have to be positive in order to get a physiologically relevant mathematical model.
The model described by ( 1) is nonnegative because it is kinetic [25], [15, pp. 153-154].This is a desirable property since the physiological problem is also nonnegative, and we will also utilize the positivity of the system in the qualitative analysis in the upcoming sections.
We use a simple control law  which is called P type control in control engineering.This produces an injection rate that is proportional to the measured tumor volume.
Suppose that k > 0 and all the model parameters are positive.Then, if x 1 increases, u increases, which results in an increase in x 3 .Increasing x 3 results in increasing Hill function x 3 ∕(ED 50 + x 3 ), which results in a decrease in x 1 .Thus, (2) provides a negative feedback for k > 0, since the tumor volume decreases as the control input increases.As a result, we will constrain k to be positive due to control considerations.

QUALITATIVE ANALYSIS OF THE CLOSED-LOOP TUMOR GROWTH MODEL
We carry out qualitative analysis of the system using computer algebra [23].Although the resulting formulas are still complex, we make some physiological considerations to simplify the calculations.For the practical problem, we can make the following assumptions: 1. a − n > 0. If this condition holds, then the tumor grows without treatment.This is the practically relevant case that should be controlled with chemotherapy.2. a − n − b < 0. If this condition holds, the chemotherapy is able to decrease tumor volume.We need this condition to be true; otherwise, the application of the therapy is not effective.
We use realistic parameter values based on measurement data acquired from mice experiments carried out in 2020 and 2021 by the Drug Resistance Research Group of the Research Center of Natural Sciences, Hungary, in collaboration with the Physiological Controls Research Center of Óbuda University.
We perform the analysis of the system using two sets of parameter values given in Tables 2 and 3.In these two sets, the pharmacokinetic and pharmacodynamic parameters have different values, showing that they can span even different orders of magnitudes.This is especially true for the pharmacokinetic parameter k 2 .The identification of the pharmacokinetic parameters is challenging from the noisy measurements.Moreover, the sensitivity of the model to the parameters ED 50 and k 2 is low; thus, relatively large differences in these parameters cause small differences in the behavior.
In Sections 3,4, and 5.1, we present the method in detail using the parameters given in Table 2.Then, we repeat the same procedure and provide the results in Section 5.2 for the parameter values given in Table 3.We note that the parameter values are given in a format of rational numbers in our initial calculations.This choice has the advantage that they can be treated symbolically, thus avoiding potential numerical errors.In later stages, however, we used high-precision decimals to enhance computational efficiency.
First, we calculate the singular points of the system.The trivial singular point is not interesting for the physiological problem since, in this case, there is no tumor.The nontrivial singular point is which is positive due to the positivity of the parameters and the assumptions made at the beginning of the section.These expressions depend on k, that is, we can control the value of the singular point with the choice of the control gain.The Jacobian at the singular point A is with ( The characteristic polynomial of the Jacobian at the singular point A is An interesting observation from a control theoretic point of view is that the characteristic polynomial does not depend on the feedback gain k.This implies that the control algorithm does not affect the qualitative properties of the singular point; these properties only depend on the parameters of the controlled system.This shows that if we can only use the measured tumor volume for the control algorithm, then we depend highly on the model parameters, which is problematic from the safety point of view since some changes in the model parameters can result in undesirable dynamic behavior, as we will show in Section 5. Now, we fix the values of w, ED 50 , c, k 1 , and k 2 as given in Table 2.With this setting, the characteristic polynomial is + 1,250, 000, 000, 000bs 4 ) . ( We find pure imaginary eigenvalues of the Jacobian using the method from [26] by looking for the parameter values that result in a characteristic polynomial which can be represented in the form with W > 0, which ensures that the polynomial in the second factor has pure imaginary roots.We look for parameter values for which the polynomial ( 8) is equal to the characteristic polynomial (7).Without loss of generality, we assume b 2 = 1.We eliminate b 0 and b 1 and get the following system for the coefficients a, b, n, and W: Setting the measured value a = 3131∕12,500 and solving Equations ( 9) and ( 10) together with the conditions using the CylindricalDecomposition command of Mathematica, we obtain suitable interval for b, that is, we get that 0.434813 ≤ b < 0.442529 or b > 0.434813.Using this result, we set the value of b and calculate the value of n accordingly.We find that the characteristic polynomial has pure imaginary roots for the parameter values where the exact value of n is a root of a fourth-degree polynomial and these are close to the identified values of the parameters in Table 2.The eigenvalues of the Jacobian with this choice of the parameters are thus, there are two pure imaginary eigenvalues.

STABILITY OF THE EQUILIBRIUM STATE AND THE HOPF BIFURCATION
As shown above for the values of parameters given in (12) and Table 2 Jacobian (4) evaluated at the singular point A has a pair of pure imaginary eigenvalues, and two other eigenvalues are negative.Thus, the system has a two-dimensional center manifold passing through A, and the singularity at the center manifold can be a stable focus, an unstable focus, or a center.In this section, we will show that, in our case, the singularity is a stable focus, and therefore, the system admits a supercritical Hopf bifurcation.A usual way to determine the type of a Hopf bifurcation is the computation of a Poincaré-Lyapunov normal form or a normal form on the center manifold, see, for example, [27,28], and it involves rather laborious computations.However, since according to the Pliss reduction principle [29], the stability of the system at the singular point A is the same as the stability of the system reduced to the center manifold, to show that the system has a stable focus at the center manifold, it is sufficient to show that A is asymptotically stable point of (1).To this end, we can use the Lyapunov functions method as described below.
According to the first Lyapunov method [22], an analytic system . where , and all eigenvalues of A ′ have negative real parts can be transformed to the normal form on the invariant surface .
where P 1 , P 2 , and Y ′ are formal series, Y ′ ( 1 ,  2 , 0) ≡ 0, P 2 = P1 , and  2 = ȳ1 .Separating the real and imaginary parts of P 1 , we write the function P 1 as If G ≡ 0, then the origin of ( 14) is stable but not asymptotically stable.When with g 0 ≠ 0, the origin is unstable if g 0 > 0, and asymptotically stable if g 0 < 0. In this case, there is a positive-defined quadratic form Q( ′ ) such that for the quadratic form it holds that for the Lie derivative of U, we have see [22, p. 87].This observation justifies the following four steps approach for the construction of a Lyapunov function for system (1) in a neighborhood of the singular point A.
the shifted system is .
Next, we choose the parameters in such a way that the Jacobian at the origin has a pair of pure imaginary eigenvalues.This is already done: The Jacobian at the origin is the same as the Jacobian at A given in ( 4) and ( 5), and if a, w, ED 50 , c, k 1 , k 2 have the values given in Table 2, and b, n are given by (12), then this condition holds.We set k = 1∕100.2. Since the right-hand sides are not polynomials, we use a Taylor series expansion of the transformed system (17) about the origin up to degree 4. Using the above parameter settings, the Taylor series expansion with approximate values is .
x 4 = 166.46x 3 − 16.738x 4 (18) For the numerical calculations, we applied the SetPrecision command of Mathematica with precision 100.It can also be verified symbolically that the Jacobian of system (18) is the same as the Jacobian J given in ( 4) and ( 5) with eigenvalues a 1 = −183.979,a 2 = −0.108065, a 3 = i, a 4 = −i where  = 0.0485016.3.As a next step, we transform system (18) in such a way that its linear part has real Jordan normal form as follows.The structure of this system is .
x = Jx +  (x) where x = (x 1 , x 2 , x 3 , x 4 ) ⊤ and  contains the nonlinear terms.If for the matrix S the product S −1 JS has real Jordan normal form, then with the substitution x = Sy where y = ( 1 ,  2 ,  3 ,  4 ) ⊤ we have that .y = S −1 .
x = S −1 JSy + S −1  (Sy).Let us denote the eigenvectors corresponding to the eigenvalues a 1 , a 2 , a 3 , a 4 by v 1 , v 2 , v 3 , v 4 , respectively.Then, v 4 is a complex conjugate of v 3 .Let S be the matrix whose columns are v 1 , v 2 , Re(v 3 ), and Im(v 3 ), respectively, then We verified numerically that Applying the transformation x = Sy for system (18) and denoting the new variables with x 1 , x 2 , x 3 , x 4 instead of  1 ,  2 ,  3 ,  4 , we get that the system has the form .
where the functions h 1 , h 2 , h 3 , h 4 contain the nonlinear terms.4.Then, we look for a polynomial of the form such that the quadratic part of Φ, is positively defined, and for the Lie derivative of the system, we have that The coefficient g 0 is called a Lyapunov coefficient or focus quantity.Comparing the coefficients of the corresponding quadratic terms on both sides of (A1), we get that It can be seen that Φ 2 will be positively defined if g 0 < 0 and C > 0. Next, we compare the corresponding third-and fourth-degree terms.With the SetPrecision command of Mathematica, we get the numerical value g 0 = −0.00878212C.It can be seen that C can be chosen arbitrarily, so with the choice C = 1, the quadratic part of Φ, is positively defined.Since g 0 is negative, it means that the origin of the transformed system and, therefore, the nontrivial singular point of the original system is locally asymptotically stable.
From the above analysis, we conclude that a supercritical Hopf bifurcation occurs at A if we perturb system (1) in such a way that Jacobian (4) has two complex conjugate eigenvalues with positive real parts.

EVOLUTION OF SYSTEM STATES UNDER PERTURBATIONS
In this section, we present numerical simulations confirming the theoretical analysis performed above and discuss the physiological meaning of the system behavior.In Section 5.1, we consider two perturbations of the system with fixed parameter values in Table 2 studied in the previous section.In Section 5.2, we provide the results of the same analysis with parameter values in Table 3.

Results with parameter values in Table 2
We choose n as a perturbation parameter and perturb its bifurcation value n 0 ≈ 0.132497 given in (12) with a positive and negative value.With the perturbation, the condition 0 < a − n < b also holds, and the value of n remains realistic.In both cases, the feedback gain is k = 1∕100.
1.If the perturbation is positive such that n = n 0 + 1∕50 ≈ 0.152497, then the singular point of the system is Since each eigenvalue has a negative real part, then the singular point is locally asymptotically stable; that is, all trajectories converge locally to the singular point.It means that the volume of the tumor tends to a constant value.
No bifurcation occurs in this case.If we increase the value of the perturbation, then the magnitude of the negative real part of the complex eigenvalues will be larger and thus the convergence of the trajectories to the singular point will be faster.2. If the perturbation is negative such that n = n 0 − 1∕50 ≈ 0.112497, then the singular point of the system is .7854, 0.08845, 0.8796) and the eigenvalues of the Jacobian are In Section 4, we showed that the Lyapunov coefficient g 0 is negative, which implies that the singular point is locally asymptotically stable.Since the Jacobian has a pair of complex eigenvalues with positive real parts, then the singular point becomes locally unstable, which means that within the region of attraction described above, on a two-dimensional central manifold, the trajectories spiral away from the singular point.From these two facts, it follows that on the central manifold, a stable limit cycle appears around the singular point due to a supercritical Hopf bifurcation.
Since the limit cycle is stable, then the trajectories approach it both from inside and from outside.In order to represent the limit cycle, we choose two initial states from the outer part and inner part of the cycle and plot the state variables as a function of the time.In these two cases, the singular point is shifted with the vectors (i) x(0) = x * + (0, 10, 10,10) ⊤ ; (ii) x(0) = x * + (0, 0, 0, 0.5) ⊤ .The results for case (i) are shown in Figure 1.The trajectory goes inward, approaching the limit cycle from outside in such a way that the amplitudes are decreasing.Figure 2 shows the results for case (ii).Here, the trajectory goes outward, approaching the limit cycle from inside in such a way that the amplitudes are increasing.The period of the oscillations, as well as the amplitudes of the trajectories close to the limit cycle, can be estimated from the figures.
In both cases, the volume of the tumor tends to the same stabilized oscillatory state, which means that it cannot be totally eliminated by the chemotherapy.The fact whether the amplitudes decrease or increase during the oscillation depends solely on the initial state of the system, that is, this behavior is independent of the control gain.
The size of the limit cycle depends on the magnitude of the perturbation.If we increase the magnitude of the negative perturbation, then the real part of the complex eigenvalues will increase; therefore, the limit cycle will be larger, and it can be observed that the trajectories will converge to it faster.
The limit cycle can also be plotted in three-dimensional phase portraits, this is shown in Figure 3 where we have chosen the state variables x 1 , x 2 , and x 3 .We remark that the system with the original parameter values given in Table 2 also contains a limit cycle, and the behavior of the system is similar to those in Figures 1-3.In order to prove this, it is enough to change one suitably chosen parameter to make the eigenvalues of the Jacobian purely imaginary.If the bifurcation parameter is reset to the original FIGURE 3 Three-dimensional phase portraits of the living tumor volume (x 1 ), the dead tumor volume (x 2 ), and the drug level in the blood (x 3 ) with k = 1∕100.In case (A), the initial condition is x(0) = x * + (0, 10, 10,10) ⊤ , and the time interval is t = 1800 days.The trajectory is going inward, approaching the limit cycle from outside.In case (B), the initial condition is x(0) = x * + (0, 0, 0, 0.5) ⊤ , and the time interval is t = 2500 days.The trajectory is going outward, approaching the limit cycle from inside.[Colour figure can be viewed at wileyonlinelibrary.com]FIGURE 4 Case (A) is a three-dimensional phase portrait of the living tumor volume (x 1 ), the dead tumor volume (x 2 ), and the drug level in the blood (x 3 ).Case (B) shows the time evolution of the living tumor volume (x 1 ).In both cases k = 1∕100, the initial condition is x(0) = x * + (0, 0, 1, 1) ⊤ , and the time interval is t = 1200 days.The trajectory is going inward, converging to the singular point.[Colour figure can be viewed at wileyonlinelibrary.com] value, then the difference between the original value and the bifurcation value corresponds to a perturbation.In (12), we have chosen b = 3∕2 close to the original value, but of course, the measured value b = 3771∕2500 would also be suitable here.In that case, the bifurcation value is n 0 ≈ 0.132558, corresponding to the perturbation −0.045915.

Results with parameter values in Table 3
The main difference between the parameter settings given in Tables 2 and 3 is that while in the former case, the system contains a stable limit cycle, in the latter case, each eigenvalue of the Jacobian has a negative real part which means that the singular point is locally asymptotically stable and thus there is no stable limit cycle around it.
Following the procedure in Section 3, we can achieve by changing the value of a suitable parameter that the Jacobian has purely imaginary eigenvalues.If we fix the values of a, w, ED 50 , c, k 1 , k 2 as given in Table 3, then we get the condition FIGURE 5 Time evolution of the living tumor volume (x 1 ) with k = 1∕100 and time interval t = 1300 days.In case (A), the initial condition is x(0) = x * + (0, 1, 1, 1) ⊤ .The trajectory is going inward, approaching the limit cycle from outside.In case (B), the initial condition is x(0) = x * + (0, 0, 0, 0.2) ⊤ .The trajectory is going outward, approaching the limit cycle from inside.[Colour figure can be viewed at wileyonlinelibrary.com] b > 1.26957.We set the measured value b = 22,229∕10,000, and with this, the bifurcation value for n is n 0 ≈ 0.0726605.If n = n 0 , then the eigenvalues of the Jacobian are s 1 = −1609.92,s 2 = −0.267993,s 3 = +0.118351i,s 4 = −0.118351i.Applying the Lyapunov method in Section 4, we get that with this parameter setting g 0 = −1.06023< 0.Then, we investigate the perturbation of n 0 both with a positive and with a negative value.In both cases, the feedback gain is k = 1∕100.
1.If we set the measured value n = 0.13881 as given in Table 3, then this corresponds to the positive perturbation n − n 0 ≈ 0.0661495.In this case, the eigenvalues of the Jacobian are Since each eigenvalue has a negative real part, then the singular point (x * 1 , x * 2 , x * 3 , x * 4 ) = (0.6444, 3.4805, 0.1816, 0.02793) is locally asymptotically stable.We plotted the phase portrait in three dimensions, choosing the state variables x 1 , x 2 , x 3 and also plotted x 1 as a function of the time starting from the initial state x(0) = x * + (0, 0, 1, 1) ⊤ , this is shown in Figure 4.If we increase the magnitude of the positive perturbation, then the convergence to the singular point will be faster.2. Next, we choose a negative perturbation such that n = n 0 − 5∕100 ≈ 0.0226605.Then, 0 < a − n < b also holds, and this value of n still falls within a range that is realistic from the viewpoint of the experiment.The singular point of the system is (x * 1 , x * 2 , x * 3 , x * 4 ) = (1.0084,5.4467, 0.2842, 0.0437), and the eigenvalues of the Jacobian are Since the complex eigenvalues of the Jacobian have positive real part and g 0 < 0, then a stable limit cycle appears around the singular point due to a supercritical Hopf bifurcation.We plotted x 1 as a function of the time starting from the initial states x(0) = x * + (0, 1, 1, 1) ⊤ and x(0) = x * + (0, 0, 0, 0.2) ⊤ ; this is shown in Figure 5.
To sum up the results of this example, we found that if we set the measured values of the parameters a, b, w, ED 50 , c, k 1 , k 2 , and n > n 0 , then the singular point is locally asymptotically stable while if n < n 0 , then a stable limit cycle appears around the singular point.
Finally, we remark that following the above procedure, it may be possible that we find a limit cycle with a parameter setting where the values of some parameters are not realistic.

CONCLUSION
We performed qualitative analysis of a tumor growth control model.We found that the system has one nontrivial singular point in the first orthant, whose coordinates depend on the feedback gain.However, the eigenvalues of the Jacobian at the singular point are independent of it, which means that the feedback gain has no impact on the stability of the singular point, and thus, it may not influence the qualitative properties of the system.
We showed that there exist realistic parameter values for which the system admits limit cycle oscillations.The qualitative method consists of the following main steps.First, setting the bifurcation value of a suitably chosen parameter, we achieved that the Jacobian has a pair of pure imaginary eigenvalues.Next, using a Lyapunov method, we showed that the first Lyapunov coefficient is negative, which means that the singular point is locally attracting.Then, perturbing the bifurcation parameter around the bifurcation value, we showed that the system can have two kinds of qualitative properties: 1.If the complex eigenvalues of the Jacobian have a negative real part, then the singular point remains locally asymptotically stable, and the trajectories of the system converge to it.In this case, the tumor volume tends to a constant value.2. If the complex eigenvalues of the Jacobian have positive real part then a stable limit cycle appears around the singular point due to a supercritical Hopf bifurcation.In this case, the tumor volume tends to a stabilized oscillatory state.The amplitudes of the oscillations either increase or decrease, which depends only on the initial state of the system.
We also represented the time evolution of the state variables.Plotting the solutions, the period of the oscillations, and the amplitudes close to the limit cycle can be estimated from the figure.Based on this, it can be decided whether the behavior of the tumor is acceptable from the medical point of view.It is favorable if the size of the tumor decreases and remains below a certain bound.The figures also show if the therapy is less efficient, this is the case when the peak value of the oscillations increases and exceeds a certain value that is no longer tolerable.
The theoretical results have two practical applications.
1.The analysis shows that a simple mathematical model, considering only a few fundamental physiological phenomena, can reproduce an oscillatory behavior, demonstrating that it is not necessary to create a complex model to describe tumor growth.Moreover, it shows medical experts that these few phenomena can lead to a repeating sequence of tumor relapse and remission.2. Controlled drug delivery systems are emerging technologies that are expected to be widely used in clinical practice in the near future.These devices allow continuous administration of the drug, like the control problem dealt with here.Our analysis showed that using a common control method, that is, injecting a dose proportional to the tumor volume, can result in a series of tumor relapses and remission.

TABLE 1
The notations, names, and dimensions of the parameters of the tumor growth model.

TABLE 2
Set #1 of the parameters based on mice experiments.

TABLE 3
Set #2 of the parameters based on mice experiments.