A fractional-order mathematical model for lung cancer incorporating integrated therapeutic approaches

This paper addresses the dynamics of lung cancer by employing a fractional-order mathematical model that investigates the combined therapy of surgery and immunotherapy. The significance of this study lies in its exploration of the effects of surgery and immunotherapy on tumor growth rate and the immune response to cancer cells. To optimize the treatment dosage based on tumor response, a feedback control system is designed using control theory, and Pontryagin’s Maximum Principle is utilized to derive the necessary conditions for optimality. The results reveal that the reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(R_0)$$\end{document}(R0) is 2.6, indicating that a lung cancer cell would generate 2.6 new cancer cells during its lifetime. The reproduction coefficient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(R_c)$$\end{document}(Rc) is 0.22, signifying that cancer cells divide at a rate that is 0.22 times that of normal cells. The simulations demonstrate that the combined therapy approach yields significantly improved patient outcomes compared to either treatment alone. Furthermore, the analysis highlights the sensitivity of the steady-state solution to variations in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_5$$\end{document}k5 (the rate of division of cancer stem cells) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_{13}$$\end{document}k13 (the rate of differentiation of cancer stem cells into progenitor cells). This research offers clinicians a valuable tool for developing personalized treatment plans for lung cancer patients, incorporating individual patient factors and tumor characteristics. The novelty of this work lies in its integration of surgery, immunotherapy, and control theory, extending beyond previous efforts in the literature.

Lung cancer is a prevalent cancer worldwide and a leading cause of cancer-related mortality 1 . The disease is characterized by genetic mutations that result in uncontrolled cell growth and division, and its development is influenced by various factors such as smoking, environmental pollutants, and genetic predisposition 2 . However, the complexity of the disease has made the development of effective treatments challenging despite advancements in cancer research 3 .
Combining PD-L1 monoclonal antibody treatment with surgery has emerged as a promising approach to treating lung cancer 4 . PD-L1 monoclonal antibodies target the PD-1/PD-L1 checkpoint pathway, which regulates the immune response 5 . These drugs block the interaction between PD-L1 on cancer cells and PD-1 on T cells, enabling the immune system to recognize and attack cancer cells 6 . Although immunotherapy is a promising treatment, prolonged administration could lead to immune-related adverse events (irAEs), including endocrine adverse events (eAEs) 7,8 . On the other hand, surgery involves the removal of the tumor and surrounding tissue and is often curative for early-stage lung cancer, but may not be effective for advanced cases 9 . Combining these therapies may enhance treatment efficacy while reducing adverse effects 10 .
Mathematical modeling has been used to study a variety of disease dynamics [11][12][13][14] and has become increasingly important in cancer research to understand cancer growth dynamics and the interactions between cancer and immune cells in the tumor microenvironment 15 . Fractional calculus, which deals with derivatives and integrals of non-integer order, has emerged as a powerful tool for modeling complex systems [16][17][18] . Fractional order differential equations have been applied in various fields such as biomedical engineering, economics, and control theory [19][20][21] .
In recent years, fractional calculus has been increasingly used in modeling cancer growth. Several mathematical models have been developed to describe the growth and spread of cancer 22 . The Caputo derivative, a common tool in fractional calculus, has been utilized to describe cancer growth dynamics and treatment response 23 . For instance, a fractional order model was used to examine the effects of chemotherapy on cancer cell growth, considering the interactions between cancer cells, chemotherapy drugs, and immune cells [24][25][26] . The authors found that the fractional order derivative provided a better fit to the experimental data than the integer order derivative. Another study utilized a fractional order model to describe breast cancer growth under the influence of immune cells, considering the interactions between cancer cells, immune cells, and chemotherapy Preliminaries Definition 1 Caputo derivative 33 The Caputo derivative of order α ∈ (0, 1] of a sufficiently differentiable function f(t) is defined as follows: where Ŵ is the gamma function.

Model formation
A fractional order system of differential equations that captures the interactions between epithelial cells, oncogenes, tumor suppressor genes, immune cells, blood vessels, and growth factors in lung cancer could take the following form: x(t) = f (x(t), u(t), t), L(x(t), u(t), t)dt, www.nature.com/scientificreports/ where: N(t) represents the number of cancer cells in the lung tissue at time t; P(t) represents the number of cancer cells that have spread to other parts of the body at time t; I(t) represents the number of immune cells in the lung tissue at time t 40 .
The model incorporates the effects of oncogenes and tumor suppressor genes through the growth rate of cancer cells, , and the carrying capacity, K as shown in Fig. 1 and equation (1). It also includes the effects of immune cells through the parameter β 1 and the growth and spread of cancer cells through the parameters µ , γ , and δ . The role of blood vessels in delivering nutrients to cancer cells and promoting metastasis is captured through the parameter β 3 . Finally, the effects of growth factors are captured through the parameters φ 1 , φ 2 , and φ 3 . Initial conditions for N(t), P(t), and I(t) must also be specified. The model can be used to simulate the growth and spread of lung cancer and investigate the effects of different treatments and interventions.

Theorem 1 The solution to model (1) exists and is unique.
Proof To prove the existence and uniqueness of solutions to the fractional order model (1), we will use the Banach fixed point theorem. Let X be the space of continuous functions on [0, T] equipped with the supremum norm | · | ∞ , where T is the time horizon. We will define an operator F : X → X and show that it is a contraction mapping on a suitable subset of X, which will imply the existence and uniqueness of a fixed point.
We define the operator F : X → X by the following equations: where R = max{N 1 0 , N 2 0 } , M = max{I 1 0 , I 2 0 , P 1 0 , P 2 0 }. We first show that F maps B R into itself. Let f ∈ B R , then we have: where we have used the fact that |f | ∞ ≤ R for all f ∈ B R . Therefore, F maps B R into itself.
We can continue the inequality by using the properties of integrals and the fact that N, I, P are bounded by K: where L = C + D = max{L 1 , L 2 , L 3 } is the Lipschitz constant. Therefore, and so F is a Lipschitz continuous function on B R with Lipschitz constant L. Now we show that F is a contraction mapping on B R , we need to show that there exists a constant 0 < A < 1 such that Since F is Lipschitz continuous on B R with constant L, using the Mean Value Theorem for integrals, we can write: where α is chosen such that (t−s) α−1 Ŵ(α) is integrable over [0, t], and M is a constant that depends on t and α. Since |Ff | ≤ RM , where, R = max{N 1 0 , N 2 0 } and M = max{I 1 0 , I 2 0 , P 1 0 , P 2 0 } , and k 1 (t) = (t−s) α−1 Ŵ(α) , k 2 (t) = (t−s) α−1 Ŵ(α) , and k 3 (t) = (t−s) α−1 Ŵ(α) are kernels of convolution integrals, using the triangle inequality and the above estimate, we get: where d = |f 1 − f 2 | ∞ . Thus, we obtain: www.nature.com/scientificreports/ from assumption. Therefore, we have shown that 0 < LM < 1 . Since A = LM then F is a contraction mapping on B R with 0 < A < 1 . This means that there exists a unique fixed point f * ∈ B R such that F(f * ) = f * . This fixed point is the solution to the system of the integral equations.
Hence, system (1) has a unique solution.

Theorem 2 System (1) is asymptotically stable at the positive equilibrium
To perform an asymptotic analysis of the fractional-order cancer model, we need to analyze the behavior of the system as time goes to infinity. In this case, we will assume that t → ∞. We will assume that the system is in a steady state when t → ∞ . Therefore, we will set all the time derivatives in the system to zero, and solve for the steady-state values of the variables N(t), I(t), and P(t).
From the first equation of system (1), we have: Since we assume that the system is not in a trivial state, we must have N(t) = 0 . Therefore, we can simplify the equation to: From the second equation of system (1), we have: From the third equation of system (1), we have: Again, since we assume that the system is not in a trivial state, we must have P(t) = 0 . Therefore, we can simplify the equation to: We can now solve this system of equations for the steady-state values of N(t), I(t), and P(t). From the equation γ N(t) − δ − β 3 I(t) = 0 , we have: Substituting this into the equation I(t)(φ 3 + β 2 P(t)) − φ 2 N(t) 2 = φ 1 I 0 , we get: Multiplying both sides by γ N(t) − δ , we get: This will lead to a cubic equation in N(t), which can be solved to obtain the steady-state value of N(t). However, since this is a complicated equation, we will not attempt to solve it explicitly.
To analyze the stability of the steady-state, we need to calculate the Jacobian matrix of the system evaluated at the steady-state values. The Jacobian matrix is given by: Evaluating the partial derivatives and substituting the steady-state values, we get: where N , Ī , and P are the steady-state values of N(t), I(t), and P(t), respectively.
To determine the stability of the steady-state, we need to calculate the eigenvalues of the Jacobian matrix. The eigenvalues can be either real or complex, and their signs determine the stability of the steady-state. If all eigenvalues have negative real parts.
To calculate the eigenvalues of the Jacobian matrix, we can we can use the characteristic equation method, which involves finding the roots of the characteristic polynomial of the matrix, given by: where is the eigenvalue and I is the identity matrix. Solving this equation for gives the eigenvalues of J.
For the current system, the characteristic equation is: where T 1 , T 2 , and T 3 are the coefficients of the characteristic polynomial. The coefficients are given by: To find the eigenvalues, we need to solve the characteristic equation, which is a cubic equation in : To show that all eigenvalues have negative real parts, we need to use the Routh-Hurwitz stability criterion. According to the criterion, all eigenvalues have negative real parts if and only if all the leading principal minors of the Routh array are positive 41 . The Routh array is constructed as follows: The first two rows of the Routh array are obtained directly from the coefficients of the characteristic equation.
The other rows are obtained by calculating the determinants of 2 × 2 submatrices of the array.  www.nature.com/scientificreports/ The conditions for the leading principal minors to be positive are: These conditions are satisfied, since we have shown that M 1 > 0 , M 2 > 0 , and M 3 > 0 . Therefore, all eigenvalues have negative real parts, which implies the system is stable.
Global stability.

Theorem 3 System (1) is globally asymptotically stable at the positive equilibrium.
Proof To prove the global stability of the system, we need to show that there exists a unique positive equilibrium point and that it is globally asymptotically stable. First, we find the equilibrium points by setting all the derivatives of (1) to zero: From the third equation of (2), we can solve for P as: Substituting this expression for P into the second equation of (2), we get: Simplifying this expression, we get a quadratic equation in I: This quadratic equation has two solutions, but we are only interested in the positive solution, which is given by: Now we substitute the expressions for I and P into the first equation of (2) and solve for N: This gives us the unique positive equilibrium point (N , I , P * ).
Next, we need to show that the equilibrium point is globally asymptotically stable. We use the Lyapunov function V (N, I, P) = N 2 + I 2 + P 2 . Taking the time derivative of V along the trajectories of the system, we get: Substituting the expressions for the equilibrium point (N * , I * , P * ) , and simplifying, we get: (2) www.nature.com/scientificreports/ Now we show that the equilibrium point is globally asymptotically stable by showing that V is decreasing along all trajectories of the system. From the expression above, we can see that if N = N * , I = I * , or P = P * , then dV dt is negative definite. If N = N * , then we have: Since I * , P * > 0 and φ 1 > φ 3 , γ N * > δ , and β 2 P * > 0 and β 3 I * > 0 , we have dV dt < 0 for all I = I * and P = P * . Similarly, if I = I * , then we have: Since N * , P * > 0 and > µP * + β 1 I * K , γ N * > δ , and β 3 I * > 0 , we have dV dt < 0 for all N = N * and P = P * . Finally, if P = P * , then we have: Since N * , I * > 0 and > µP * + β 1 I * K , φ 1 > φ 3 , and β 2 P * > 0 , we have dV dt < 0 for all N = N * and I = I * . Therefore, we have shown that dV dt < 0 along all trajectories of the system that do not coincide with the equilibrium point (N * , I * , P * ) . Since V is a Lyapunov function and is decreasing along all trajectories, the equilibrium point is globally asymptotically stable.
Reproduction number and coefficient. The reproduction number and reproduction coefficient are used to measure the potential of an infectious disease to spread in a population. In this case, we can adapt these concepts to measure the potential of cancer to grow and spread in the body.
To compute the reproduction number R 0 and the reproduction coefficient R c of system (1), we need to first determine the disease-free equilibrium point E 0 = (N 0 , I 0 , P 0 ) . This is the point at which all populations are in their uninfected or baseline state, i.e., N 0 = K , I 0 = 0 , and P 0 = 0.
To compute the Jacobian matrix of the system evaluated at E 0 , we differentiate each equation with respect to N, I, and P: Evaluating these partial derivatives at E 0 yields: where J(E 0 ) is the Jacobian matrix evaluated at the disease-free equilibrium E 0 .
Next, we need to compute the eigenvalues of J(E 0 ) , which will give us information about the stability of E 0 . If all eigenvalues have negative real parts, then E 0 is stable and the disease (i.e., cancer) will not grow or spread. If at least one eigenvalue has a positive real part, then E 0 is unstable and the disease has the potential to grow and spread.
The eigenvalues of J(E 0 ) are given by the roots of the characteristic polynomial: To compute the reproduction number and reproduction coefficient, we need to determine the dominant eigenvalue of J(E 0 ) . If this eigenvalue is real and positive, then the disease has the potential to grow and spread.
Let 1 be the dominant eigenvalue of J(E 0 ) . Then the reproduction number is given by: and the reproduction coefficient is given by:

Sensitivity analysis.
To perform sensitivity analysis, we will use the concept of the normalized sensitivity coefficient 42,43 , defined as: where S i is the sensitivity coefficient for the ith parameter, and ln denotes the natural logarithm. The sensitivity coefficient measures the proportional change in N(∞) resulting from a proportional change in k i . A sensitivity coefficient of S i = 1 means that a 1 To calculate the sensitivity coefficients, we first need to find the steady-state solution of the system, which satisfies: Substituting the expressions for each of the variables, we obtain the following equations: Solving these equations for C * , P * , S * , E * , N * , and I * , we get: www.nature.com/scientificreports/ Next, we will calculate the sensitivity coefficients S i using the expression: where ∂N(∞) ∂k i is the partial derivative of N(∞) with respect to k i . Next, we will calculate the sensitivity of the solution with respect to the parameters using the formula: where θ is a parameter and ǫ is a small perturbation. We will use a finite difference approximation to compute this derivative, with ǫ = 0.01θ.
First, we will define the parameters and the initial conditions: Next, we will compute the sensitivities of the solutions with respect to the parameters: Taking the derivative of N(t) with respect to k 11 , we have: Using the equations for Ṅ (t) and N(t) from earlier, we can simplify this expression as: Using the expressions for Ṅ (t) and k 12 from earlier, we get: Simplifying this expression, we get: Taking the derivative of I(t) with respect to k 14 , we have: Using the equations for İ (t) and I(t) from earlier, we can simplify. The sensitivity of the solution with respect to parameter k i is given by: . 11 . 11 . www.nature.com/scientificreports/ where δ ij is the Kronecker delta, which is equal to 1 if i = j and 0 otherwise. This expression shows how the solution u(t) changes as we perturb the value of each parameter k i . The terms inside the square brackets represent the contributions of each parameter to the overall sensitivity of the solution.

Control theory
To apply control theory to the fractional order lung cancer model, we will first need to define our control objective. In this case, we will aim to design a control strategy that minimizes the population of infected cancer cells (I(t)) while minimizing the use of resources (P(t)).
To accomplish this objective, we will use a Proportional-Integral-Derivative (PID) controller to generate control signals that will modulate the population of resources and infected cancer cells in the model. Specifically, we will use the error between the desired and actual population of infected cancer cells as feedback to compute a control signal that will regulate the populations of both infected cancer cells and resources.
The control objective can be expressed as follows: where K p , K i , and K d are the proportional, integral, and derivative gains of the controller, respectively, and e(t) = I d (t) − I(t) is the error signal.
To simplify the notation, we will use the following definitions: www.nature.com/scientificreports/ where f 1 , f 2 , and f 3 are the right-hand side functions of the first-order fractional differential equations, and u 1 and u 2 are the control signals for the infected cancer cells and resources, respectively. We set u 2 to zero since we only want to control the population of infected cancer cells and resources. Using the definitions above, we can express the control problem as: Next, we will apply the Pontryagin's Maximum Principle to derive the necessary conditions for optimality. Let x(t) = [N(t), I(t), P(t)] be the state vector of the system, and u(t) = [u 1 (t), u 2 (t)] be the control signal vector. The Hamiltonian of the system is defined as: is the costate vector, and L(x(t), u(t)) is the cost function to be minimized.
The necessary conditions for optimality are: Stationarity condition: Transversality condition: State dynamics: Costate dynamics: Boundary conditions: where t 0 and t f are the initial and final time, respectively, and x 0 is the initial state. Note that the stationarity condition gives the optimal control signal u 1 (t) in terms of the costate variable 2 (t) , which can be computed by solving the state and costate equations numerically. Once 2 (t) is computed, the optimal control signal can be obtained by substituting it back into the stationarity condition.
The optimal control problem can be solved by formulating the above equations as an initial value problem and using numerical methods to solve it. The solution will provide the optimal trajectory for the system variables N(t), I(t), and P(t), as well as the optimal control signal u 1 (t).

Combined therapy optimization.
To incorporate PDL1 monoclonal antibody immunotherapy and surgery as a control for the lung cancer model, we can modify the third equation to include the effects of the treatment. Specifically, we can assume that the treatment reduces the population of cancer cells and tumorpromoting immune cells, which can be modeled by adding a term proportional to the product of the treatment Minimize: www.nature.com/scientificreports/ efficacy and the populations of cancer cells and immune cells. We can also assume that the treatment has a time delay before it becomes effective. The modified system of equations is: where T(t) is the population of tumor cells that have been surgically removed, θ(t − τ ) is the Heaviside step function with a time delay τ , ǫ 1 , ǫ 2 , ǫ 3 , and ǫ 4 are constants representing the efficacy of the treatment and the surgery, and ǫ 1 and ǫ 2 correspond to the reduction in the populations of cancer cells and immune cells, respectively, due to the treatment. The term θ(t − τ )ǫ 3 P(t) represents the influx of tumor cells back into the body from the surgical site, which can be assumed to occur with a delay of τ after surgery. The term −ǫ 4 T(t) represents the removal of the surgically removed tumor cells from the body.
This modified system of equations can be used to study the effects of the combination of PDL1 monoclonal antibody immunotherapy and surgery on the growth of lung cancer cells. The efficacy of the treatment and the surgery can be varied to explore their impact on the growth of the cancer cells and the outcome of the treatment.
To compute the coefficient of the modified system of equations, we need to identify the terms that contain the variables of interest and the parameters that affect their interaction. The variables of interest are N(t), I(t), P(t), and T(t). The parameters that affect their interaction are µ , β 1 , φ 2 , φ 3 , β 2 , γ , δ , β 3 , ǫ 1 , ǫ 2 , ǫ 3 , and ǫ 4 . The product of the treatment efficacy and the populations of cancer cells and immune cells is given by We differentiate the first equation with respect to N(t) to obtain: To compute the coefficient of I(t), we need to differentiate the second equation with respect to I(t). We obtain: Therefore, the coefficient of To compute the coefficient of P(t), we need to differentiate the second, third, and fourth equations with respect to P(t). We obtain: To compute the coefficient of T(t), we need to differentiate the fourth equation with respect to T(t). We obtain: To develop an objective function to improve treatment, dosage and reduce cancer growth in this model, we can consider the following goals: Minimize the population of cancer cells, N(t). Minimize the population of tumor-promoting immune cells, I(t). Minimize the population of tumor cells that have not been surgically removed, P(t).
Maximize the population of surgically removed tumor cells, T(t).
To achieve these goals, we can formulate the following objective function: www.nature.com/scientificreports/ where w 1 , w 2 , w 3 , and w 4 are weighting factors that represent the relative importance of each goal. The integral is taken over the time horizon [0, t f ] , where t f is the final time of the simulation. The first term, w 1 N(t) , represents the goal of minimizing the population of cancer cells. The second term, w 2 I(t) , represents the goal of minimizing the population of tumor-promoting immune cells. The third term, w 3 P(t) , represents the goal of minimizing the population of tumor cells that have not been surgically removed. The fourth term, −w 4 T(t) , represents the goal of maximizing the population of surgically removed tumor cells.
We can use this objective function to guide the optimization of treatment dosage and timing, and to explore the impact of different treatment and surgery parameters on the outcome of the treatment. By adjusting the weighting factors, we can place more or less emphasis on each goal, depending on the specific priorities of the patient and the medical team.
Applying Pontryagin's Maximum Principle in (6), we introduce the Hamiltonian function as follows: where ψ i (t) are the costate variables, which represent the rate of change of the cost function with respect to the state variables.
Using the Hamiltonian, we can write the necessary conditions for optimality as follows: Now, we need to solve the system of differential equations that we obtained from the Pontryagin's Maximum Principle. We start with the first equation: where N 0 is the initial condition of N.

Next, we solve for T(t):
Integrating both sides of the equation with respect to time, we obtain: where T 0 is the initial condition of T. Next, we solve for I(t): To find I(t), we first notice that the term θ(t − τ ) is equal to 1 when t ≥ τ and equal to 0 when t < τ . Therefore, we split the solution into two cases: Case 1: t < τ Case 2: t ≥ τ For Case 1, the equation simplifies to: To solve for I(t) in this case, we need to use the boundary condition I(τ ) = I 0 . We can write the solution as: For Case 2, the equation simplifies to: To solve for I(t) in this case, we need to use the boundary condition I(τ ) = I 0 . We can write the solution as: Next, we solve for V(t), which represents the vascular endothelial growth factor (VEGF) : Integrating both sides of the equation with respect to time, we obtain: where V 0 is the initial condition of V. Finally, we solve for W(t), which represents the wound healing cells: Integrating both sides of the equation with respect to time, we obtain: where W 0 is the initial condition of W.

Feedback control.
To design a feedback control system that adjusts the treatment dosage based on the tumor response, we can use the modified system of equations: (21) I(t) = I(τ )e − t τ (β 1 ψ 2 +φ 3 ψ 2 +β 2 ψ 5 +β 3 ψ 7 +ǫ 2 ψ 9 )ds . www.nature.com/scientificreports/ where ǫ 1 , ǫ 2 , ǫ 3 , and ǫ 4 are the constants representing the efficacy of the treatment and the surgery, and θ(t − τ ) is the Heaviside step function with a time delay τ representing the delay before the treatment becomes effective.
To design the feedback control system, we need to define a control variable that can be adjusted based on the tumor response to the treatment. One possible control variable is the population of cancer cells N(t), which we can use to adjust the dosage of the treatment. We can define a feedback control system that adjusts the treatment dosage based on the difference between the target population of cancer cells N target and the actual population of cancer cells N(t).
The feedback control system can be implemented as follows: 1. Set an initial dosage for the treatment. The PID controller can be defined as follows: The proportional term (P) adjusts the treatment dosage based on the difference between the target and the actual population of cancer cells. The larger the difference, the larger the adjustment in the dosage. The integral term (I) accumulates the error over time and adjusts the treatment dosage to reduce the accumulated error. The derivative term (D) adjusts the treatment dosage based on the rate of change of the error, which can help to prevent overshooting the target population of cancer cells. Repeat steps 2-4 until the actual population of cancer cells is within a specified tolerance of the target population.
The target population of cancer cells N target can be set based on the stage and severity of cancer, as well as the treatment goals and desired outcomes. The target population of cancer cells may also need to be adjusted based on the patient's response to the treatment.
To obtain the time evolution of the populations of cancer cells, immune cells, and surgically removed tumor cells for the current treatment dosage, we would set the fractional-order to be 1 and solve the resulting system of equations:

where N(t) is the population of cancer cells, I(t) is the population of immune cells, P(t) is the population of surgically removed tumor cells, and T(t) is the population of tumor cells that were not removed by surgery.
The explicit solution for the system of differential equations is: where, The explicit solution provides a mathematical expression for each compartment as a function of time, given the initial conditions and the values of the model parameters. We can calculate the difference between the target population of cancer cells N target and the actual population of cancer cells N(t) at the end of the treatment period. If the actual population of cancer cells is within a specified tolerance of the target population, we can stop the treatment. Otherwise, we can adjust the dosage of the treatment based on the difference between the target and the actual population of cancer cells using a proportional-integral-derivative (PID) controller.
The PID controller can be defined as follows: where u(t) is the treatment dosage, e(t) is the error between the target population of cancer cells N target and the actual population of cancer cells N(t), K p , K i , and K d are the proportional, integral, and derivative gains, respectively. The proportional term (P) adjusts the treatment dosage based on the difference between the target and actual population of cancer cells. The larger the difference, the larger the adjustment in the dosage. The integral term (I) accumulates the error over time and adjusts the treatment dosage to reduce the accumulated error. The derivative term (D) adjusts the treatment dosage based on the rate of change of the error, which can help to prevent overshooting the target population of cancer cells.
The PID controller can be used to adjust cancer treatment dosage until the actual population of cancer cells is within a specified tolerance of the target population. The target population can be set based on the cancer stage, severity, treatment goals, and desired outcomes. The PID controller uses the time evolution of populations to calculate the error and adjust the treatment dosage. The process is repeated until the actual population of cancer cells is within the specified tolerance.
To use the PID controller to adjust the dosage of the treatment, we need to calculate the error e(t), which is the difference between the target population of cancer cells N target and the actual population of cancer cells N(t) at time t. The error can be expressed as: We can then use the PID controller to compute the treatment dosage u(t) as: where K p , K i , and K d are the proportional, integral, and derivative gains, respectively.
To implement the PID controller for the cancer treatment model, we can follow the following steps: Step 1: Set the initial values for the population of cancer cells N 0 and surgically removed tumor cells P 0 . Set the target population of cancer cells N target and the specified tolerance.
Step 2: Set the initial values for the integral and derivative terms of the PID controller, I 0 and D 0 , respectively.
Step 3: At each time step t, solve the system of differential equations using the explicit solution and obtain the actual population of cancer cells N(t) and surgically removed tumor cells P(t).
Step 4: Calculate the error e(t) = N target − N(t).
Update the integral and derivative terms of the PID controller: Integral term: Step 5: Compute the treatment dosage using the PID controller: Proportional term: P(t) = K p e(t).
If the absolute value of the error |e(t)| is less than the specified tolerance, stop the treatment. Otherwise, continue to the next time step.

Results and discussion
Analytical results show that the system has a unique positive equilibrium point, and is globally stable. The dynamics of the lung cancer model without treatment can be seen in Fig. 2. The simulation provides insights into the population dynamics of cancer cells (N), cancer cells that have spread (P), and immune cells (I) in the context of lung cancer. The different fractional order values ( α ) of 0.5, 0.7, and 0.9 demonstrate varying effects on these populations. A notable observation is that higher fractional order values lead to faster convergence of the model, indicating improved accuracy and efficiency in describing the cancer dynamics. As ( α ) increases, the growth rate of normal cells (N) decreases, suggesting a stronger inhibitory effect on cell proliferation. This implies that higher fractional order values enable the model to capture the control mechanisms governing cancer cell growth within lung tissue more effectively. Furthermore, the plot reveals that increasing ( α ) results in a increase in the population of cancer cells that have spread (P), indicating an increased ability of cancer cells to form secondary tumors in other parts of the body. The faster convergence associated with higher fractional order values  www.nature.com/scientificreports/ contributes to a more accurate representation of the suppression of metastasis, emphasizing the importance of incorporating a strong fractional derivative component in modeling the immune response against spreading cancer cells. In contrast, the immune cell population (I) exhibits fluctuating patterns, with oscillations becoming more pronounced as ( α ) increases. This indicates that higher fractional order values lead to a more dynamic immune response, resulting in periodic fluctuations in the number of immune cells. The faster convergence of the model with higher fractional order values allows for a more precise depiction of the intricate interplay between immune cells and cancer cells, highlighting the complex dynamics of the immune system's interaction with lung cancer. The sensitivity coefficients for each parameter were calculated using the concept of the normalized sensitivity coefficient. The analysis indicated that the most sensitive parameter is the carrying capacity (K), followed by the growth rate of cancer cells ( ) and the rate at which cancer cells spread from the lung tissue to other parts   www.nature.com/scientificreports/ of the body ( µ ), as shown in Table 2 and Fig. 7. This suggests that interventions aimed at reducing the carrying capacity or inhibiting the growth or spread of cancer cells could be effective in reducing the number of cancer cells in the lung tissue and preventing their spread to other parts of the body. The analysis revealed that the parameter with the greatest impact on the steady-state solution of the system is the rate of tumor cell growth, k 1 . Increasing k 1 by 1 percent would result in a 78.5 percent increase in the steady-state tumor cell concentration while decreasing the effectiveness of treatment by reducing k 2 by 1 percent would lead to a 130.7 percent increase in tumor cell concentration. The analysis also showed that the steady-state solution is most sensitive to changes in the rate of division of cancer stem cells ( k 5 ) and the rate of differentiation of cancer stem cells into progenitor cells ( k 13 ). Interventions that target cancer stem cells, such as therapies that inhibit stem cell division or promote their differentiation, may have a significant impact on the overall tumor cell population. The effects of immune cells, blood vessels, and growth factors were found to be less sensitive, indicating that these factors may play a less significant role in the growth and spread of lung cancer. However, this does not mean that they are not important, as they may still have important roles in specific stages or types of lung cancer. Overall, the sensitivity analysis provides valuable insights into the key parameters that drive the growth and spread of lung cancer and can inform the development of targeted interventions. Furthermore, results show that the reproduction number R 0 represents the average number of cancer cells that are produced by a single cancer cell over the course of its lifetime. A value of 2.6 for lung cancer would indicate that, on average, each cancer cell would produce 2.6 (at least 2) new cancer cells during its lifetime as shown in Fig. 3. This can contribute to the rapid growth and spread of cancer within the body. The reproduction coefficient R C represents the rate at which cancer cells reproduce or divide, relative to the rate at which normal cells in the body reproduce or divide. A value of 0.22 for lung cancer would indicate that cancer cells are dividing at a rate that is 0.22 times the rate at which normal cells divide. This could suggest that cancer cells are growing more slowly than normal cells in the body, but it is important to note that cancer growth is a complex and multifactorial process that involves many different factors beyond just cell division.
The result shows that the optimization models improve treatment, and dosage and reduce cancer growth. The objective function minimizes the population of cancer cells, tumor-promoting immune cells, and tumor cells that have not been surgically removed, and maximizes the population of surgically removed tumor cells, formulated using the weighting factors as seen in Fig. 4. The lung cancer optimization combined therapy model describes a modified system of equations that incorporates the effects of PDL1 monoclonal antibody immunotherapy and surgery as a control for lung cancer. The modified system of equations includes a term proportional to the product of treatment efficacy and the populations of cancer cells and immune cells to model the reduction of these populations due to the treatment. Additionally, the treatment has a time delay before it becomes effective.
The PID controller consists of three terms -proportional, integral, and derivative gains -each of which adjusts the dosage based on different factors. The proportional term adjusts the dosage based on the error between the target and the actual population of cancer cells, with larger errors resulting in larger adjustments. The integral term accumulates the error over time and adjusts the dosage to correct for consistent differences between the target and actual population. The derivative term adjusts the dosage based on the rate of change of the error, preventing overshooting of the target population. The magnitude of the adjustments for each term is proportional to the gain constant for that term -K p for the proportional term, K i for the integral term, and K d for the derivative term. By using the time evolution of the populations and the PID controller to adjust the dosage, the actual population of cancer cells can be brought within a specified tolerance of the target population, with the feedback control shown in Fig. 5. The target population of cancer cells can be set based on the stage and severity www.nature.com/scientificreports/ of the cancer, treatment goals, and desired outcomes, and may need to be adjusted based on patient response and changes in the cancer over time.
In Fig. 6, the W(t) plot demonstrates an initial surge in wound healing cells followed by a gradual decline, indicating the body's acute response to tissue injury and subsequent resolution of the healing process. This suggests that the lung tissue's healing capacity is dynamic, with an initial boost in wound-healing cells that gradually subsides as healing progresses. On the other hand, the v(t) plot reveals a sustained elevation of vascular endothelial growth factor, indicating a continuous need for angiogenesis to support tissue repair. This implies that the lung tissue requires a prolonged vascular supply to facilitate proper healing.
The employed approach has demonstrated high efficiency. The use of fractional order derivatives in this lung cancer model provides relevance by incorporating memory and non-local dependencies into the system. The fractional order ( α ) influences the rate of change and interactions between different cell populations, representing the memory effects in cellular processes. By including fractional derivatives, the model captures the long-term behavior and complex interactions involved in lung cancer dynamics, contributing to a more accurate representation of the disease's progression and potential therapeutic interventions. The novelty of this research lies in the development of a comprehensive fractional-order mathematical model for lung cancer that incorporates integrated therapeutic approaches. This approach allows for a more precise understanding of the intricate interactions among cancer cells, immune cells, and other constituents of the tumor microenvironment. Additionally, by incorporating PD-L1 monoclonal antibody treatment and surgery as controls in the model, this study can explore the potential benefits of combining these treatments and investigate their impact on tumor dynamics. This approach stands out from others because it offers a more comprehensive and accurate representation of lung cancer dynamics and treatment response. It also holds promise for guiding the development of novel treatment strategies for lung cancer and enhancing patient outcomes.