To numerically solve systems (1), (2), and (10), we consider the initial value problem in (1):$$\begin{aligned} {_0^c}D_t^{\alpha }{\textbf{y}}(t)={\textbf{f}}(t, {\textbf{y}}(t)), \quad {\textbf{y}}(0)={\textbf{y}}_0. \end{aligned}$$Employing the Riemann–Liouville integral operator in Definition 2.2, we get that:$$\begin{aligned} {\textbf{y}}(t)-{\textbf{y}}_0=\frac{1}{\Gamma (\alpha )}\int _0^t (t-s)^{\alpha -1}\, {\textbf{f}}(s, {\textbf{y}}(s))\, ds. \end{aligned}$$
(24)
After substituting \(t=t_{n+1}\) into Eq. (24) and subtracting two obtained equations, we can write:$$\begin{aligned} {\textbf{y}}(t_{n+1})={\textbf{y}}_0+\frac{1}{\Gamma (\alpha )} \sum \limits _{m=0}^n \int _{t_m}^{t_{m+1}} (t_{n+1}-s)^{\alpha -1}\, {\textbf{f}}(s, {\textbf{y}}(s))\, ds, \end{aligned}$$
(25)
where \(t_j=jh, j=0, 1, \ldots , N\) and \(h=T_f/N\) is the step size. We now approximate the function \({\textbf{f}}(s, {\textbf{y}}(s))\) on the interval \([t_m, t_{m+1}]\) using the two-step Lagrange polynomial interpolation:$$\begin{aligned} {\textbf{f}}(s, {\textbf{y}}(s))&\approx \frac{s-t_{m+1}}{t_m-t_{m+1}}{\textbf{f}}(t_m, {\textbf{y}}_m)+ \frac{s-t_{m}}{t_{m+1}-t_{m}}{\textbf{f}}(t_{m+1}, {\textbf{y}}_{m+1})\\&=-\frac{s-t_{m+1}}{h}{\textbf{f}}(t_m, {\textbf{y}}_m)+ \frac{s-t_{m}}{h}{\textbf{f}}(t_{m+1}, {\textbf{y}}_{m+1}), \end{aligned}$$
(26)
where \({\textbf{y}}_{k}={\textbf{y}}(t_k)\). Using (25) and (26), we have$$\begin{aligned} {\textbf{y}}_{n+1}&={\textbf{y}}_0 + \frac{1}{h\, \Gamma (\alpha )} \bigg \{ \sum \limits _{m=0}^{n} \int _{t_m}^{t_{m+1}} (t_{n+1}-s)^{\alpha -1}\, (s-t_m)\, {\textbf{f}}(t_{m+1}, {\textbf{y}}_{m+1})\, ds \\&\quad -\sum \limits _{m=0}^{n} \int _{t_m}^{t_{m+1}} (t_{n+1}-s)^{\alpha -1}\, (s-t_{m+1})\, {\textbf{f}}(t_{m}, {\textbf{y}}_{m})\, ds \bigg \}, \quad n=0, 1, \ldots , N.\end{aligned}$$
(27)
Using integration by parts, (27) is converted into the following formula:$$\begin{aligned} {\textbf{y}}_{n+1}&={\textbf{y}}_0 + \frac{h^{\alpha }}{\Gamma (\alpha +2)} \sum \limits _{m=0}^{n} \bigg \{ [(n-m+1)^{\alpha +1}-(n-m)^{\alpha }(n-m+\alpha +1)]\, {\textbf{f}}(t_{m+1}, {\textbf{y}}_{m+1}^p) \\&\quad -[(n-m+1)^{\alpha }(n-m-\alpha )-(n-m)^{\alpha +1}]\, {\textbf{f}}(t_{m}, {\textbf{y}}_{m})\bigg \}, \quad n=0, 1, \ldots , N. \end{aligned}$$
(28)
Due to appearing \({\textbf{y}}_{m+1}\) in the right side of (28), this formula is an implicit formula and values of \({\textbf{y}}_{m+1}\) should be predicted (as \({\textbf{y}}_{m+1}^p\)). Thus, formula (28) will be a corrector formula. In formula (25), we use the rectangle rule for the integral part and obtain the following predictor formula:$$\begin{aligned} {\textbf{y}}_{n+1}^p ={\textbf{y}}_{0}+\frac{h^{\alpha }}{\Gamma (\alpha +1)} \sum \limits _{m=0}^{n} B_{\frac{n(n+1)}{2}+m+1}\,{\textbf{f}}(t_{m}, {\textbf{y}}_{m}), \quad n=0, 1, \ldots , N, \end{aligned}$$
(29)
where,$$\begin{aligned} B_{\frac{n(n+1)}{2}+m+1}=(n-m+1)^{\alpha }-(n-m)^{\alpha }, \quad n=0, 1, \ldots , N, \,\, m=0, 1, \ldots , n. \end{aligned}$$Therefore, the numerical formula for system (1) is as follows:The predictor formula:$$\begin{aligned} N_{n+1}^p&=N_0+\frac{h^{\alpha }}{\Gamma (\alpha +1)} \sum \limits _{m=0}^n B_{\frac{n(n+1)}{2}+m+1}\,\bigg \{ \lambda ^{\alpha }N_m\left( 1 – \frac{N_m}{K^{\alpha }}\right) – \mu ^{\alpha } N_m P_m – \beta ^{\alpha }_1N_m I_m – \zeta ^{\alpha }_1 N_m M_m\bigg \}, \\ I_{n+1}^p&=I_0+ \frac{h^{\alpha }}{\Gamma (\alpha +1)} \sum \limits _{m=0}^n B_{\frac{n(n+1)}{2}+m+1}\,\bigg \{ \phi ^{\alpha }_1I_0 + \phi ^{\alpha }_2N_m^2 – \phi ^{\alpha }_3I_m – \beta ^{\alpha }_2I_m P_m + \zeta ^{\alpha }_2 I_m M_m – \zeta ^{\alpha }_3 I_m R_m\bigg \},\\ P_{n+1}^p&=P_0+ \frac{h^{\alpha }}{\Gamma (\alpha +1)} \sum \limits _{m=0}^n B_{\frac{n(n+1)}{2}+m+1}\,\bigg \{ \gamma ^{\alpha } N_m P_m – \delta ^{\alpha } P_m – \beta ^{\alpha }_3I_m P_m + \zeta ^{\alpha }_4 P_m M_m\bigg \},\\ M_{n+1}^p&=M_0 +\frac{h^{\alpha }}{\Gamma (\alpha +1)} \sum \limits _{m=0}^n B_{\frac{n(n+1)}{2}+m+1}\,\bigg \{\zeta ^{\alpha }_5 N_m M_m – \zeta ^{\alpha }_6 I_m M_m – \zeta ^{\alpha }_7 P_m M_m\bigg \},\\ R_{n+1}^p&=R_0+\frac{h^{\alpha }}{\Gamma (\alpha +1)} \sum \limits _{m=0}^n B_{\frac{n(n+1)}{2}+m+1}\,\bigg \{\zeta ^{\alpha }_8 I_m R_m \bigg \}. \end{aligned}$$The corrector formula:$$\begin{aligned} N_{n+1}&=N_0+\frac{h^{\alpha }}{\Gamma (\alpha +2)} \sum \limits _{m=0}^n \bigg [ ((n-m+1)^{\alpha +1}-(n-m)^{\alpha }(n-m+\alpha +1))\bigg \{ \lambda ^{\alpha }N^p_{m+1}\left( 1 – \frac{N^p_{m+1}}{K^{\alpha }}\right) \\&\quad – \mu ^{\alpha } N_{m+1}P_{m+1} ^p- \beta ^{\alpha }_1N_{m+1}^p I_{m+1}^p – \zeta ^{\alpha }_1 N_{m+1}^p M_{m+1}^p \bigg \}\\&\quad -((n-m+1)^{\alpha }(n-m-\alpha )-(n-m)^{\alpha +1})\bigg \{ \lambda ^{\alpha }N_m\left( 1 – \frac{N_m}{K^{\alpha }}\right) – \mu ^{\alpha } N_mP_m – \beta ^{\alpha }_1N_mI_m – \zeta ^{\alpha }_1 N_mM_m \bigg \} \bigg ], \\ I_{n+1}&=I_0+\frac{h^{\alpha }}{\Gamma (\alpha +2)} \sum \limits _{m=0}^n \bigg [ ((n-m+1)^{\alpha +1}-(n-m)^{\alpha }(n-m+\alpha +1))\bigg \{ \phi ^{\alpha }_1I_0 + \phi ^{\alpha }_2{N^2}_{m+1}^p – \phi ^{\alpha }_3I_{m+1}^p\\&\quad – \beta ^{\alpha }_2I_{m+1}^p P_{m+1}^p + \zeta ^{\alpha }_2 I_{m+1}^p M_{m+1}^p – \zeta ^{\alpha }_3 I_{m+1}^p R_{m+1}^p\bigg \}-((n-m+1)^{\alpha }(n-m-\alpha )-(n-m)^{\alpha +1})\bigg \{ \phi ^{\alpha }_1I_0 \\&\quad + \phi ^{\alpha }_2N_m^2 – \phi ^{\alpha }_3I_m – \beta ^{\alpha }_2I_m P_m + \zeta ^{\alpha }_2 I_m M_m – \zeta ^{\alpha }_3 I_m R_m \bigg \} \bigg ], \\ P_{n+1}&=P_0+\frac{h^{\alpha }}{\Gamma (\alpha +2)} \sum \limits _{m=0}^n \bigg [ ((n-m+1)^{\alpha +1}-(n-m)^{\alpha }(n-m+\alpha +1))\bigg \{ \gamma ^{\alpha } N_{m+1}^p P_{m+1}^p – \delta ^{\alpha } P_{m+1}^p\\&\quad – \beta ^{\alpha }_3I_{m+1}^p P_{m+1}^p + \zeta ^{\alpha }_4 P_{m+1}^p M_{m+1}^p \bigg \}-((n-m+1)^{\alpha }(n-m-\alpha )-(n-m)^{\alpha +1})\bigg \{ \gamma ^{\alpha } N_m P_m\\&\quad – \delta ^{\alpha } P(t) – \beta ^{\alpha }_3I_m P_m + \zeta ^{\alpha }_4 P_m M_m \bigg \} \bigg ], \\ M_{n+1}&=M_0+\frac{h^{\alpha }}{\Gamma (\alpha +2)} \sum \limits _{m=0}^n \bigg [ ((n-m+1)^{\alpha +1}-(n-m)^{\alpha }(n-m+\alpha +1))\bigg \{ \zeta ^{\alpha }_5 N_{m+1}^p M_{m+1}^p – \zeta ^{\alpha }_6 I_{m+1}^p M_{m+1}^p\\&\quad – \zeta ^{\alpha }_7 P_{m+1}^p M_{m+1}^p\ \bigg \}-((n-m+1)^{\alpha }(n-m-\alpha )-(n-m)^{\alpha +1})\bigg \{ \zeta ^{\alpha }_5 N_m M_m – \zeta ^{\alpha }_6 I_m M_m – \zeta ^{\alpha }_7 P_m M_m\ \bigg \} \bigg ],\\ R_{n+1}&=R_0+\frac{h^{\alpha }}{\Gamma (\alpha +2)} \sum \limits _{m=0}^n \bigg [ ((n-m+1)^{\alpha +1}-(n-m)^{\alpha }(n-m+\alpha +1))\bigg \{ \zeta ^{\alpha }_8 I_{m+1}^p R_{m+1}^p \bigg \}\\&\quad -((n-m+1)^{\alpha }(n-m-\alpha )-(n-m)^{\alpha +1})\bigg \{ \zeta ^{\alpha }_8 I_m R_m \bigg \} \bigg ]. \end{aligned}$$Optimal system (14) can be solved by above algorithm. Similarly, we can use the following method to solve system (15):$$\begin{aligned} \varvec{\lambda }_{n+1}&=\ \frac{h^{\alpha }}{\Gamma (\alpha +2)} \sum \limits _{m=0}^{n} \bigg \{ [(n-m+1)^{\alpha +1}-(n-m)^{\alpha }(n-m+\alpha +1)]\, \varvec{\Lambda }(t_{m+1}, {\textbf{y}}_{m+1}, \varvec{\lambda }_{m+1}^p) \\&\quad -[(n-m+1)^{\alpha }(n-m-\alpha )-(n-m)^{\alpha +1}]\, \varvec{\Lambda }(t_{m}, {\textbf{y}}_{m}, \varvec{\lambda }_{m})\bigg \}, \\ \varvec{\lambda }_{n+1}^p&=\frac{h^{\alpha }}{\Gamma (\alpha +1)} \sum \limits _{m=0}^{n} B_{\frac{n(n+1)}{2}+m+1}\,\varvec{\Lambda }(t_{m}, {\textbf{y}}_{m}, \varvec{\lambda }_{m}), \quad n=0, 1, \ldots , N, \end{aligned}$$where \(\varvec{\Lambda }(t, {\textbf{y}}(t), \varvec{\lambda }(t))\) is the vector of equations in the right side of system (15). After solving systems (1), (14), and (15), values of control variables can be updated by (16).For numerical simulation, the following values have been considered for parameters and initial conditions in system (1):$$\begin{aligned} \lambda&=0.3, \,\, \mu =0.01, \,\, \beta _1=0.01, \,\, \zeta _1=0.9, \,\, \phi _1=0.03, \,\, \phi _2=0.004, \,\, \phi _3=0.1, \,\, \beta _2=0.01, \,\, \zeta _2=0.458, \,\, K=10000, \,\,\\ \zeta _3&=0.8, \,\, \gamma =0.003, \,\, \delta =0.03, \,\, \beta _3=0.4, \,\, \zeta _4=0.045, \,\, \zeta _5=0.5, \,\, \zeta _6=0.6, \,\, \zeta _7=0.5, \,\, \zeta _8=0.6, \\ N_0&=1.5, \,\, I_0=0.7, \,\, P_0=0.3, \,\, M_0=0.1, \,\, R_0=0.1. \end{aligned}$$Figure 3 depicts the behaviour of the state variables in system (1) for \(\alpha =0.7, 0.8, 0.9, 0.95, 1\). Figure 4 depicts the behaviour of the state variables in system (1) for \(\alpha =0.7, 0.8, 0.9, 0.95, 1\), \(N_0=1, I_0=3, P_0=1, M_0=1, R_0=1\). The number of cancer cells, spread cancer cells, and enhanced immune cells (in millions) increases and the number of immune cells starts to decrease after increasing. In all cases, the figures plotted for diverse values of \(\alpha\) approach the figure plotted for \(\alpha =1\). In Fig. 5, actual data points are compared to predicted values obtained from the suggested algorithms and the model. A coincidence between real data and numerical values is seen over the interval [0, 5] and on [5, 15], the simulated figures have an increasing or decreasing behaviour similar to the figures of the real data. Figures of state variables in optimal system (14)-(17) are seen in Fig. 12 for \(\alpha =0.7, 0.8, 0.9, 0.95, 1\), \(N_0=1, I_0=3, P_0=1, M_0=1, R_0=1\). In order to survey the validity of the numerical results obtained from the suggested model, the absolute residual errors for the state variables in System (1) are calculated. For this purpose, all terms in the equations of System (1) are shifted to the left side and obtained numerical values are substituted into them:$$\begin{aligned} {\mathcal {R}}_1(t)&={_0^c}D_{t}^{\alpha }N(t) -\lambda ^{\alpha }N(t)\left( 1 – \frac{N(t)}{K^{\alpha }}\right) + \mu ^{\alpha } N(t)P(t) + \beta ^{\alpha }_1N(t)I(t) + \zeta ^{\alpha }_1 N(t)M(t) \approx 0,\\ {\mathcal {R}}_2(t)&={_0^c}D_{t}^{\alpha }I(t)- \phi ^{\alpha }_1I_0 – \phi ^{\alpha }_2N(t)^2 + \phi ^{\alpha }_3I(t) + \beta ^{\alpha }_2I(t)P(t) – \zeta ^{\alpha }_2 I(t)M(t) + \zeta ^{\alpha }_3 I(t)R(t)\approx 0,\\ {\mathcal {R}}_3(t)&={_0^c}D_{t}^{\alpha }P(t) -\gamma ^{\alpha } N(t)P(t) + \delta ^{\alpha } P(t) + \beta ^{\alpha }_3I(t)P(t) – \zeta ^{\alpha }_4 P(t)M(t)\approx 0,\\ {\mathcal {R}}_4(t)&= {_0^c}D_{t}^{\alpha }M(t)+\zeta ^{\alpha }_5 N(t)M(t) + \zeta ^{\alpha }_6 I(t)M(t) + \zeta ^{\alpha }_7 P(t)M(t)\approx 0,\\ {\mathcal {R}}_5(t)&={_0^c}D_{t}^{\alpha }R(t)- \zeta ^{\alpha }_8 I(t)R(t)\approx 0, \end{aligned}$$where \({\mathcal {R}}_i(t), i=1, 2, 3, 4, 5\) are residual functions. Thus, we have the following system for \(\alpha =1\):$$\begin{aligned} {\mathcal {R}}_1(t_j)&=\frac{N_{j+1}-N_j}{h} -\lambda ^{\alpha }N_j\left( 1 – \frac{N_j}{K^{\alpha }}\right) + \mu ^{\alpha } N_j P_j + \beta ^{\alpha }_1N_j I_j + \zeta ^{\alpha }_1 N_j M_j \approx 0,\\ {\mathcal {R}}_2(t_j)&=\frac{I_{j+1}-I_j}{h}- \phi ^{\alpha }_1I_0 – \phi ^{\alpha }_2N_j^2 + \phi ^{\alpha }_3I_j + \beta ^{\alpha }_2I_jP_j – \zeta ^{\alpha }_2 I_j M_j + \zeta ^{\alpha }_3 I_j R_j\approx 0,\\ {\mathcal {R}}_3(t_j)&=\frac{P_{j+1}-P_j}{h} -\gamma ^{\alpha } N_j P_j + \delta ^{\alpha } P_j + \beta ^{\alpha }_3I_j P_j – \zeta ^{\alpha }_4 P_j M_j \approx 0,\\ {\mathcal {R}}_4(t_j)&= \frac{M_{j+1}-M_j}{h}+\zeta ^{\alpha }_5 N_j M_j + \zeta ^{\alpha }_6 I_j M_j + \zeta ^{\alpha }_7 P_j M_j \approx 0,\\ {\mathcal {R}}_5(t_j)&=\frac{R_{j+1}-R_j}{h}- \zeta ^{\alpha }_8 I_j R_j\approx 0, \end{aligned}$$for \(j=0, 1, \cdots , N\). Figures of absolute residual errors are depicted in Fig. 6. As expected, the values of the residual functions are small. In other words, the obtained numerical results from the suggested algorithms are getting close to the exact data if they are available. As another criterion to measure the validity of the proposed model, the sensitivity of the model to some of its parameters is investigated. Hence, the plots of state variables can be seen in Figs. 7, 8, 9, 10 and 11 for \(\alpha =0.8\), initial conditions \(N_0=1.5, I_0=0.7, P_0=0.3, M_0=0.1, R_0=0.1\), and different values of diverse parameters. As can be seen, by increasing the values of parameters, figures of state functions are not divergent. In Figure 7 by increasing values of \(\beta _1\) and \(\mu\), the number of cancer cells are decreasing gradually. In Fig. 8, by increasing the value of \(\beta _2\), the number of immune cells sounds constant, while by increasing the value of \(\phi _3\), the number of immune cells decreases gradually. Cancer cells spread in a similar way by increasing values of \(\delta\) and \(\zeta _4\) in Fig. 9. In Fig. 10, no variation observes in the behaviour of M(t) by increasing the value of \(\zeta _6\). In Fig. 11, with the increase of the value of \(\zeta _8\), the number of enhanced cytotoxic immune cells increases. The figures of the control variables \(D_I(t)\) and \(D_T(t)\) are depicted in Fig. 13. The number of cancers and the spread of cancer cells is increasing with time (weeks).Figure 3Plots of cancer cells, immune cells, spread cancer cells, genetic mutations, and enhanced immune cells in system (1) for different values of \(\alpha\) and \(N_0=1.5, I_0=0.7, P_0=0.3, M_0=0.1, R_0=0.1\).Figure 4Plots of cancer cells, immune cells, spread cancer cells, genetic mutations, and enhanced immune cells in system (1) for different values of \(\alpha\) and \(N_0=1, I_0=3, P_0=1, M_0=1, R_0=1\).Figure 5Validation comparison plot of model (1) versus synthetic data.Similarly, to estimate values of \(D_T(t), D_I(t), u_T(t)\), and \(u_I(t)\) in system (18), we consider the following Hamiltonian function:$$\begin{aligned} {\textsf{H}}(t)=&a_1^{\alpha } N(t)+a_2^{\alpha } I(t)+a_3^{\alpha } P(t)+ a_4^{\alpha } M(t)+a_5^{\alpha } R(t)+b_1^{\alpha } D_I^2(t) +b_2^{\alpha } D_T^2(t) +c_1^{\alpha } u_I^2(t)+c_2^{\alpha } u_T^2(t)\\&+\Lambda _1(t) \bigg \{ \lambda ^{\alpha }N(t)\left( 1 – \frac{N(t)}{K^{\alpha }}\right) – \mu ^{\alpha } N(t)P(t) – \beta ^{\alpha }_1N(t)I(t) – \zeta ^{\alpha }_1 N(t)M(t) – \chi ^{\alpha }_1 N(t)\left( D_I(t) + u_I(t)\right) \bigg \} \\&+\Lambda _2(t) \bigg \{ \phi ^{\alpha }_1I_0 + \phi ^{\alpha }_2N^2(t) – \phi ^{\alpha }_3I(t) – \beta ^{\alpha }_2I(t)P(t) + \zeta ^{\alpha }_2 I(t)M(t) – \zeta ^{\alpha }_3 I(t)R(t) + \chi ^{\alpha }_2 I(t)\left( D_T(t) + u_T(t)\right) \bigg \} \\&+\Lambda _3(t) \bigg \{ \gamma ^{\alpha } N(t)P(t) – \delta ^{\alpha } P(t) – \beta ^{\alpha }_3I(t)P(t) + \zeta ^{\alpha }_4 P(t)M(t) – \chi ^{\alpha }_3 P(t)\left( D_T(t) + u_T(t)\right) \bigg \}\\&+\Lambda _4(t) \bigg \{ \zeta ^{\alpha }_5 N(t)M(t) – \zeta ^{\alpha }_6 I(t)M(t) – \zeta ^{\alpha }_7 P(t)M(t) + \chi ^{\alpha }_4 M(t)(D_T(t) + u_T(t)) + \chi ^{\alpha }_5 M(t)\left( D_I(t) + u_I(t)\right) \bigg \}\\&+ \Lambda _5(t) \bigg \{\zeta ^{\alpha }_8 I(t)R(t) + \chi ^{\alpha }_6 R(t)\left( D_I(t) + u_I(t)\right) \bigg \}, \end{aligned}$$
(30)
where \(\Lambda _i(t), i=1, \ldots , 5\) are adjoint variables. If \(D^*_T, D^*_I, u^*_T,\) and \(u^*_I\) are optimal values of control variables, then the optimal system, utilizing Hamiltonian (30), will be as follows:$$\begin{aligned} _{0}^{c}D_{t}^{\alpha }N(t)&= \lambda ^{\alpha }N(t)\left( 1 – \frac{N(t)}{K^{\alpha }}\right) – \mu ^{\alpha } N(t)P(t) – \beta ^{\alpha }_1N(t)I(t) – \zeta ^{\alpha }_1 N(t)M(t) – \chi ^{\alpha }_1 N(t)\left( D^*_I(t) + u^*_I(t)\right) , \nonumber \\ _{0}^{c}D_{t}^{\alpha }I(t)&= \phi ^{\alpha }_1I_0 + \phi ^{\alpha }_2N^2(t) – \phi ^{\alpha }_3I(t) – \beta ^{\alpha }_2I(t)P(t) + \zeta ^{\alpha }_2 I(t)M(t) – \zeta ^{\alpha }_3 I(t)R(t) + \chi ^{\alpha }_2 I(t)\left( D^*_T(t) + u^*_T(t)\right) , \nonumber \\ _{0}^{c}D_{t}^{\alpha }P(t)&= \gamma ^{\alpha } N(t)P(t) – \delta ^{\alpha } P(t) – \beta ^{\alpha }_3I(t)P(t) + \zeta ^{\alpha }_4 P(t)M(t) – \chi ^{\alpha }_3 P(t)\left( D^*_T(t) + u^*_T(t)\right) , \nonumber \\ _{0}^{c}D_{t}^{\alpha }M(t)&= \zeta ^{\alpha }_5 N(t)M(t) – \zeta ^{\alpha }_6 I(t)M(t) – \zeta ^{\alpha }_7 P(t)M(t) + \chi ^{\alpha }_4 M(t)(D^*_T(t) + u^*_T(t)) + \chi ^{\alpha }_5 M(t)\left( D^*_I(t) + u^*_I(t)\right) , \nonumber \\ _{0}^{c}D_{t}^{\alpha }R(t)&= \zeta ^{\alpha }_8 I(t)R(t) + \chi ^{\alpha }_6 R(t)\left( D^*_I(t) + u^*_I(t)\right) , \end{aligned}$$
(31)
$$\begin{aligned} \frac{d \Lambda _1(t)}{dt}&=-a_1^{\alpha }-\Lambda _1(t)\bigg [\lambda ^{\alpha } -\frac{2 \lambda ^{\alpha }}{K^{\alpha }}N(t)-\mu ^{\alpha }P(t)-\beta _1^{\alpha } I(t)-\zeta _1^{\alpha } M(t)-\chi _1^{\alpha } (D_I(t)+u_I(t)\bigg ]\nonumber \\&\quad -2 \phi _2^{\alpha } \Lambda _2(t) N(t)-\gamma ^{\alpha } \Lambda _3(t) P(t)-\zeta _5^{\alpha } \Lambda _4(t) M(t), \nonumber \\ \frac{d \Lambda _2(t)}{dt}&=-a_2^{\alpha }+\beta _1^{\alpha } \Lambda _1(t)N(t)-\Lambda _2(t)\bigg [ -\phi _3^{\alpha }-\beta _2^{\alpha } P(t)+\zeta _2^{\alpha } M(t)-\zeta _3^{\alpha } R(t)+\chi _2^{\alpha } (D_T(t)+u_T(t)) \bigg ]\nonumber \\&\quad + \beta _3^{\alpha } \Lambda _3(t) P(t) +\zeta _6^{\alpha } \Lambda _4(t) M(t) -\zeta _8^{\alpha } \Lambda _5(t) R(t),\nonumber \\ \frac{d \Lambda _3(t)}{dt}&=-a_3^{\alpha }+\mu ^{\alpha } \Lambda _1(t) N(t) +\beta _2 ^{\alpha } \Lambda _2(t) I(t)-\Lambda _3(t) \bigg [ \gamma ^{\alpha } N(t)-\delta ^{\alpha } -\beta _3^{\alpha } I(t)+\zeta _4^{\alpha } M(t)-\chi _3^{\alpha } (D_T(t)\nonumber \\&\quad +u_T(t)) \bigg ]+\zeta _7^{\alpha } \Lambda _4(t) M(t), \nonumber \\ \frac{d \Lambda _4(t)}{dt}&=-a_4^{\alpha }+\zeta _1^{\alpha } \Lambda _1(t) N(t) -\zeta _2^{\alpha } \Lambda _3(t) P(t)-\Lambda _4(t)\bigg [ \zeta _5^{\alpha } N(t)-\zeta _6^{\alpha } I(t)-\zeta _7^{\alpha } P(t) +\chi _4^{\alpha } (D_T(t)\nonumber \\&\quad +u_T(t))+\chi _5^{\alpha } (D_I(t)+u_I(t))\bigg ],\nonumber \\ \frac{d \Lambda _5(t)}{dt}&=-a_5^{\alpha }+\zeta _3^{\alpha } \Lambda _2(t) I(t)-\Lambda _5(t) \bigg [ \zeta _8^{\alpha } I(t) +\chi _6^{\alpha } (D_I(t)+u_I(t))\bigg ], \end{aligned}$$
(32)
$$\begin{aligned} D_I^*&=\min \{\max \{0, \Delta _I\}, 1\}, \quad D_T^*=\min \{\max \{0, \Delta _T\}, 1\},\nonumber \\ u_I^*&=\min \{\max \{0, \Psi _I\}, 1\}, \quad u_T^*=\min \{\max \{0, \Psi _T\}, 1\}, \end{aligned}$$
(33)
where$$\begin{aligned} \Delta _I&=\frac{\chi _1^{\alpha } \Lambda _1(t) N(t)-\chi _5^{\alpha } \Lambda _4(t) M(t) -\chi _6^{\alpha } \Lambda _5(t) R(t) }{2 b_1 ^{\alpha }},\\ \Delta _T&=-\frac{\chi _2^{\alpha } \Lambda _2(t) I(t)-\chi _3^{\alpha } \Lambda _3(t) P(t) +\chi _4^{\alpha } \Lambda _4(t) M(t) }{2 b_2 ^{\alpha }},\\ \Psi _I&=\frac{\chi _1^{\alpha } \Lambda _1(t) N(t)-\chi _5^{\alpha } \Lambda _4(t) M(t) -\chi _6^{\alpha } \Lambda _5(t) R(t) }{2 c_1 ^{\alpha }},\\ \Psi _T&=-\frac{\chi _2^{\alpha } \Lambda _2(t) I(t)-\chi _3^{\alpha } \Lambda _3(t) P(t) +\chi _4^{\alpha } \Lambda _4(t) M(t) }{2 c_2 ^{\alpha }}. \end{aligned}$$After solving problem (31)–(33) using the proposed predictor-corrector method, figures of state and control variables are depicted in Figs. 14 and 15. The number of spread cancer cells remains almost constant after a decreasing trend. The behaviour of control signals (\(u_I(t)\) and \(u_T(t)\)) after adjusting the drug dosages (\(D_I(t)\) and \(D_T(t)\)) is seen in Fig. 16.The model’s long-term effects E(t), defined by evolution equation (22), are seen in Fig. 17 for \(\alpha =0.7, 0.75, 0.85, 0.95, 1\), \(\theta _1=1, \theta _2=0.3, \theta _3=0.5\), and \(E_0=1\).
Figures of the quality of life QoL(t) introduced in (23) are seen in Fig. 18 for different values of parameters and \(\alpha\) and \(Q_0=1\).Now, by having approximate values of \(D_I, D_T\), And E, we can compute values of the direct costs (\(C_{\text {direct}}\)), indirect costs (\(C_{\text {indirect}}\)), and cost-benefit ratio (CBR) for \(h=0.01\) and \(T_f=3\) by the Trapezoidal method to evaluate the integral in (19) and (20). Values of these quantities are listed in Table 3 for different values of \(\alpha\), \(\delta _0=0.6\), and \(\eta =0.25\).Table 3 Direct costs, indirect costs, and cost-benefit ratio for different values of \(\alpha\)
A fractional-order model for optimizing combination therapy in heterogeneous lung cancer: integrating immunotherapy and targeted therapy to minimize side effects
