\documentclass[11pt]{article} \setlength{\topmargin}{-0.3in} \renewcommand{\baselinestretch}{1.5} \setlength{\topskip}{0.1in} \setlength{\textheight}{9.2in} \setlength{\oddsidemargin}{0.1in} \setlength{\evensidemargin}{0.1in} \setlength{\textwidth}{6.75in} \usepackage{amsmath} \allowdisplaybreaks \usepackage{float,amsthm} \usepackage{booktabs} \usepackage{amssymb} \usepackage{amsfonts} \usepackage[colorlinks=true,citecolor=blue,urlcolor=blue]{hyperref} \usepackage[mathlines,displaymath]{lineno} \runninglinenumbers \usepackage[numbers,comma, sort&compress]{natbib} %\usepackage{natbib} %\newcommand\Mycite[1]{% %\citeauthor{#1}~[\citeyear{#1}]} \usepackage{float} \usepackage{epsfig} \usepackage{epstopdf} \usepackage[toc]{appendix} \usepackage{placeins} \usepackage{graphicx} \newcommand{\R}{\mathbb R} \newcommand{\C}{\mathbb C} \newtheorem{corollary}{Corollary} \newcommand{\ds}{\displaystyle} \newtheorem{theorem}{Theorem}[section] \newtheorem{prop}{Proposition}[section] \newtheorem{definition}{Definition}[section] \newtheorem{remark}{Remark}[section] \newtheorem{lemma}{Lemma}[section] %\usepackage{mathtools} \numberwithin{equation}{section} \begin{document} \title{\bf A mathematical model of lung cancer incorporating drug resistance, macrophage polarization, and immune escape} \vspace{6in} \author{{\normalsize Willie K. Chidaushe$^{1} ,$ ~~Thomas Musora$^{1} ,$ ~~Steady Mushayabasa$\,^{2}$ \footnote{Correspondence email.: mhelikumi@yahoo.co.uk (M.H); steadymushaya@gmail.com (SM)}} \\ {\scriptsize \it $^{1}\,$ Department of Mathematics and Statistics, School of Natural Sciences and Mathematics,}\\ {\scriptsize \it Chinhoyi University of Technology, Chinhoyi, Zimbabwe }\\ {\scriptsize \it $^{2}\,$ Department of Mathematics \& Computational Sciences, University of Zimbabwe,} \\ {\scriptsize \it P.O. Box MP 167 Mount Pleasant, Harare, Zimbabwe}} %\vspace{0.1in}} \date{} \maketitle \begin{abstract} \noindent Lung cancer remains a leading cause of cancer-related mortality globally, presenting a major clinical challenge across both smoking and non-smoking populations. Deciphering the complex architecture of the tumor microenvironment, specifically the dynamic cross-talk between drug-resistant malignant sub-populations, effector cells, and infiltrating macrophages, is crucial for designing optimal therapeutic regimens that mitigate treatment failure. In this study, we formulate a data-driven, non-linear system of ordinary differential equations (ODEs) to rigorously quantify cellular evolutionary dynamics, M1/M2 macrophage polarization, and immune escape mechanisms during lung cancer progression. The mathematical validation of the model confirms its epidemiological and biological well-posedness, establishing that state solutions remain strictly non-negative and bounded within a defined invariant region. Using the next generation matrix operator, we analytically derive the basic reproduction number ($R_0$), a pivotal threshold governing cellular proliferation and tumor persistence. Local and global sensitivity analyses were performed to identify the neoplastic proliferation rate, macrophage phenotypic transition rates, and therapeutic efficacy parameters as the primary drivers governing $R_0$, thereby offering quantitative guidance for adaptive drug dosing strategies. Model selection and comparison form a critical component of this study. To identify the most parsimonious and accurate model, we evaluate candidate models using both the Akaike Information Criterion, small sample AIC and the Bayesian Information Criterion. \end{abstract} \noindent\textbf{Keywords:} Lung cancer; Mathematical modelling; Drug resistance, Tumour-immune interactions; Macrophage polarization; Immune escape \noindent\textbf{MSC 2020:} 92C50; 34D20; 34D23; 49J15; 92B05; 92C60. \section{Introduction} The interplay between drug resistance and macrophage polarization represents a multifaceted relationship that profoundly influences disease outcomes. When macrophage polarization is dysregulated, it can pave the way for the development and proliferation of drug-resistant pathogens, including notorious strains like Methicillin-resistant Staphylococcus aureus (MRSA) \cite{yan2025progress}. This emergence poses significant challenges for effective clinical treatment, complicating the management of infections. Therefore, gaining a deeper understanding of the mechanisms that govern macrophage polarization and its intricate role in drug resistance is crucial. This knowledge is key to developing innovative, targeted therapies that can enhance treatment efficacy and ultimately improve patient outcomes \cite{li2025mitochondrial}.\\ Macrophages are a vital component of the innate immune system, playing a crucial role in defending the body against foreign invaders. Their discovery dates back to 1882 when the renowned immunologist Elia Metchnikoff first identified these remarkable cells in the hatchlings of starfish. Using tangerine tree thistles as part of his study, Metchnikoff demonstrated how these cells actively participate in the process of phagocytosis, effectively engulfing and digesting harmful particles \cite{chen2023macrophages}. His groundbreaking research continued with observations of macrophages in Daphnia magna, commonly known as the water flea. In these tiny crustaceans, which were infested with various fungal spores, he highlighted the essential function of macrophages in identifying and eliminating threats, thereby maintaining the organism's health \cite{jeyachandran2025invertebrate}. This groundbreaking work laid the foundation for our understanding of macrophages as key players in the immune response, responsible for cleaning up foreign materials and facilitating the body’s defenses against infections \cite{chen2023macrophages, jeyachandran2025invertebrate}. \\ Macrophages are essential components of the immune system, playing a crucial role in regulating inflammation and modulating immune responses. These adaptable cells can modify their functions based on the specific characteristics of their surrounding microenvironment. This remarkable flexibility enables macrophages to respond effectively to various physiological and pathological conditions, ensuring a customized immune reaction that is vital for maintaining homeostasis and fighting infections \cite{he2026macrophage}. The concept of macrophage polarization vividly illustrates the remarkable adaptability of macrophages, a vital type of immune cell in the body. These cells can transform into different functional types, predominantly classified as M1 and M2 macrophages. M1 macrophages are known for their role in initiating and amplifying pro-inflammatory responses, acting as frontline defenders against pathogens and infections. They release a variety of inflammatory cytokines that help mobilize other immune cells to the site of injury or infection \cite{khayati2023potential}.\\ On the other hand, M2 macrophages take on a different role characterized by their involvement in anti-inflammatory responses. These cells are crucial for promoting healing and tissue repair after an injury. They secrete growth factors and cytokines that help resolve inflammation and facilitate the regeneration of damaged tissues \cite{li2024rabeprazole}. Striking a delicate balance between M1 and M2 macrophages is essential for maintaining the body’s homeostasis. An excessive dominance of M1 macrophages can lead to chronic inflammation and tissue damage, while an overabundance of M2 macrophages may impair the body’s ability to effectively combat infections. Thus, the interplay between these two macrophage types is vital for managing inflammation and ensuring overall health \cite{khayati2023potential}. Maintaining a proper balance between these two types is crucial for keeping the body in a state of homeostasis and effectively managing inflammation. When this balance is disrupted, it can lead to various health issues, including chronic inflammatory diseases, metabolic disorders, and cancer \cite{khayati2023potential, li2024rabeprazole}.\\ In a healthy lung, macrophages play a crucial role in maintaining tissue integrity and safeguarding against abnormal cells, effectively acting as the body's first line of defense. However, the dynamics change dramatically in the presence of lung cancer cells . These malignant cells have the ability to manipulate incoming monocytes, transforming them into Tumor-Associated Macrophages (TAMs) \cite{solinas2009tumor}. This deliberate re-education causes the macrophages to undergo a significant shift in their function; they transition from an anti-tumor phenotype, characterized by their M1-like properties, which actively combat tumors, to a pro-tumorigenic and immunosuppressive phenotype akin to M2-like behavior \cite{dunsmore2024timing}. As a result of this transformation, the TAMs become a formidable protective barrier for the tumor, effectively shielding it from various treatments. This alteration in macrophage function enables the cancer cells to develop resistance to diverse therapeutic approaches, complicating treatment efforts and contributing to the persistence and progression of the disease across all major classes of lung cancer treatment \cite{solinas2009tumor}.\\ The primary causes of drug resistance in cancer therapy can be traced back to specific genetic mutations that lead to an overproduction of drug transporters, which are proteins that facilitate the removal of therapeutic agents from cells. This over expression hampers the efficacy of chemotherapy and targeted therapies, rendering them less effective in eliminating cancerous cells \cite{ashrafi2022current}. Moreover, tumors are composed of a diverse array of cells, each exhibiting a unique combination of genetic, epigenetic, and phenotypic characteristics. This inherent heterogeneity means that different tumor cells can respond in various ways to the same treatment, which can lead to the survival of certain clones that develop resistance \cite{jamal2017tracking}. Additionally, the tumor microenvironment comprised of surrounding cells, blood vessels, and signaling molecules can undergo adaptive changes in response to treatment. These alterations can create a survival advantage for cancer cells, allowing them to thrive and proliferate even in the face of medical intervention. Together, these factors not only complicate treatment efforts but also underscore the dynamic and evolving nature of cancer as it adapts to thwart therapeutic strategies \cite{dunsmore2024timing}. \subsection{Immune Escape} Immune escape in lung cancer is a significant factor that limits the effectiveness of immunotherapy. This phenomenon involves several mechanisms, including metabolic reprogramming, overexpression of immune checkpoint molecules, and abnormalities in the presentation of antigens \cite{wang2025lung}. These mechanisms allow tumors to evade detection and response from the immune system, resulting in treatment resistance and recurrence of the disease. Recent studies have identified new therapeutic targets and promising clinical trials aimed at addressing these immune escape mechanisms. The goal is to develop personalized treatment strategies that can successfully overcome the challenges posed by immune escape in lung cancer \cite{li2026advancing}.\\ Immune escape, also known as immune evasion or antigenic escape, occurs when the immune system fails to recognize or eliminate a pathogen or cancerous cell. This phenomenon is crucial in the context of infections, cancer progression, and vaccine resistance, as it enables the invader to persist despite the body’s immune defenses or therapeutic interventions \cite{wang2025lung}. This process can occur in various ways related to both genetic and environmental factors. Such mechanisms include homologous recombination, as well as manipulation and resistance of the host's immune responses \cite{hanada2014genetic, wang2025lung}. Tumors from various types of cancer utilize strategies to evade destruction by the host immune system. Immune escape has been observed in solid tumors, including small cell and non-small cell lung cancers. These escape strategies often involve tumor-induced alterations in the tumor microenvironment (TME), making it more challenging for the immune system to recognize and respond to tumor cells effectively \cite{tufail2025immune}.\\ Furthermore, immune escape mechanisms empower malignant cells to systematically bypass host immunosurveillance, thereby accelerating subclonal selection, survival, and uncontrolled proliferation. This phenotypic adaptation encompasses a spectrum of distinct biological pathways, including somatic genetic alterations that yield neoantigenic heterogeneity, the creation of an immunosuppressive tumor microenvironment via inhibitory cytokine secretion, and the down-regulation of major histocompatibility complex class I (MHC-I) molecules to avoid direct recognition by cytotoxic T lymphocytes \cite{wang2025lung}. These dynamic adaptations severely compromise the therapeutic efficacy of conventional chemotherapies, targeted inhibitors, and cancer vaccines, presenting a major barrier to long-term disease remission and therapeutic durability \cite{hanada2014genetic}. Consequently, elucidating the precise quantitative and biochemical frameworks governing immune evasion is paramount for advancing oncology, particularly in optimizing immune checkpoint blockade regimens, engineering CAR-T cell platforms, and designing novel combination immunotherapies capable of overriding resistance mechanisms \cite{tufail2025immune}. In this research initiative, we will undertake a detailed exploration of the complex relationship between drug resistant strains of various pathogens and the polarization of macrophages key immune cells that are essential for initiating and regulating immune responses . Macrophages, which can adopt different functional states depending on the signals they receive from their environment, play a significant role in determining the outcome of infections and the effectiveness of immunotherapies \cite{chen2023macrophages}. Polarization of these cells can enhance or suppress the immune response to pathogens, and understanding how drug-resistant strains influence this polarization is critical \cite{yan2025progress}. In addition to investigating macrophage behavior, we will also analyze the various mechanisms used by pathogens to evade immune detection, including genetic mutations, antigenic variation, and the secretion of immunomodulatory factors that disrupt normal immune function \cite{tufail2025immune}.\\ To systematically investigate these non-linear population dynamics, we formulate a compartmental system of ordinary differential equations (ODEs) that serves as the deterministic framework for simulating the temporal interactions between drug-resistant malignant cell lineages and infiltrating macrophage populations \cite{cartelle2019computational}. By conducting rigorous qualitative and quantitative analyses of the model’s trajectory solutions, we elucidate the fundamental mechanistic pathways governing subclonal selection, M1/M2 macrophage polarization, and microenvironmental immune evasion. Specifically, this mathematical architecture enables the identification of critical biological tipping points and parameter sensitivities, establishing a rational baseline for designing multi-target therapeutic interventions. Ultimately, this modeling approach aims to optimize adaptive drug dosing protocols, mitigate treatment failure driven by acquired resistance, and provide actionable insights for advancing personalized therapeutic strategies in clinical oncology \cite{oliveira2024overview}. \section{Methods} Mathematical oncology provides a rigorous quantitative framework for modeling the spatio-temporal dynamics and evolutionary complexity of lung cancer. By translating complex cellular, biophysical, and clinical processes into systems of differential equations, researchers can systematically evaluate tumor progression and therapeutic response. Central to these frameworks are continuum population kinetics and deterministic ordinary differential equations (ODEs), which enable global compartmental modeling of distinct cellular subpopulations including sensitive and resistant tumor strains, effector immune cells, and polarized macrophages \cite{yan2025progress}. In this section, we utilize an ODE-based modeling architecture augmented with phenomenological growth functions designed to capture the physical and microstructural constraints imposed by the lung parenchyma. These spatial and anatomical limitations govern cellular proliferation rates, localized mechanical stress, and spatial spread. Our mathematical framework explicitly incorporates microvascular transport dynamics, which govern the transvascular exchange and intra-tumoral convective diffusive transport of essential metabolic substrates such as glucose and dissolved oxygen alongside systemic chemotherapeutic agents into the poorly perfused tumor core. Incorporating microvascular delivery kinetics is biophysically essential, as chaotic tumor angiogenesis engenders profound spatial heterogeneity, elevated interstitial fluid pressure, and localized hypoxia \cite{oliveira2024overview}. This hypoxic state serves as a primary biological driver of the epithelial mesenchymal transition, genomic instability, and secondary drug resistance. To accurately capture the non-linear, rate limited cellular and biochemical interactions within the heterogeneous tumor microenvironment, we integrate non-linear enzyme kinetics and saturating functional responses, specifically drawing upon Holling Type II and Michaelis-Menten conceptual formulations. These saturating response functions reflect fundamental physiological capacity constraints across three critical axes. Modeling the non-linear cytotoxic capacity of effector immune cells such as cytotoxic T-lymphocytes and natural killer cells—prevents non-physical, unbounded cell elimination rates at high tumor burdens, effectively capturing physiological effector cell exhaustion and target cell crowding \cite{tufail2025immune}. \\ Accounting for the saturating concentration thresholds of microenvironmental signaling networks including interferon-gamma, interleukin-4, interleukin-10, and transforming growth factor beta that drive the phenotypic transition of anti-tumor M1 macrophages into pro-tumorigenic, immunosuppressive M2 phenotypes. Quantifying the capacity-limited binding of therapeutic agents to cell-surface receptors, such as PD-1/PD-L1 immune checkpoint axes and targeted tyrosine kinase inhibitors, thereby accounting for drug-target saturation kinetics and concentration-dependent clearance pathways \cite{chen2023macrophages}. By coupling parenchymal spatial limitations and microvascular perfusion kinetics with non-linear saturating biochemical responses within a unified differential framework, this approach offers deep mechanistic resolution into the co-evolutionary dynamics of lung cancer progression. Crucially, this biophysical realism establishes a rigorous foundation for formulating optimal control strategies to evaluate adaptive drug dosing schedules that minimize systemic toxicity while actively preventing the subclonal selection of drug-resistant tumor strains. \subsection{Model derivation} We formulate a mechanistic mathematical model describing the interactions between lung tumor cells, the immune system, macrophage phenotypes, and chemotherapy. The model extends the tumor--immune interaction framework of Eftimie et al.~\cite{Eftimie2021} by explicitly incorporating drug-sensitive and drug-resistant tumor cell populations together with macrophage polarization. The model aims to investigate the emergence of chemotherapy resistance and the role of the tumor microenvironment in regulating tumor progression. The total tumor population is divided into drug-sensitive and drug-resistant tumor cells. Drug-sensitive tumor cells proliferate, compete for limited resources, are eliminated by chemotherapy and immune cells, and may acquire resistance due to treatment-induced selective pressure and M2 macrophage-mediated immunosuppression. Drug-resistant tumor cells possess reduced sensitivity to chemotherapy and continue proliferating despite treatment. The immune compartment consists of cytotoxic CD8$^{+}$ T cells together with two macrophage phenotypes. M1 macrophages exhibit anti-tumor activity by directly killing tumor cells, whereas M2 macrophages promote tumor growth and suppress cytotoxic immune responses. Chemotherapy concentration is represented using a one-compartment pharmacokinetic model. The dynamics of the model are governed by the following system of nonlinear ordinary differential equations. \begin{subequations} \label{eq:model} \begin{align} \frac{dT_S}{dt}=&r_ST_S\left(1-\frac{T_S+T_R}{K_T}\right)+\alpha M_2T_S-\delta_1M_1T_S-\kappa_1ET_S-\psi C T_S-\mu(M_2,C)T_S+\eta T_R,\label{eq:TS}\\[2ex] \frac{dT_R}{dt}=&r_RT_R\left(1-\frac{T_S+T_R}{K_T}\right)+\alpha M_2T_R-\delta_2M_1T_R-\kappa_2ET_R+\mu(M_2,C)T_S-\eta T_R-\varepsilon\psi C T_R,\label{eq:TR} \\[2ex] \frac{dE}{dt}=&\Lambda_E+\frac{\rho_0(T_S+T_R)}{a_0+T_S+T_R}E-\mu_EE-\phi M_2E-\theta E^2,\label{eq:E}\\[2ex] \frac{dM_1}{dt}=&p_{M1}M_1\left(1-\frac{M_1+M_2}{K_M}\right)+\frac{\rho_1EM_2}{a_1+E}-\mu_{M_1}M_1-\delta(T_S+T_R)M_1,\label{eq:M1}\\[2ex] \frac{dM_2}{dt}=&p_{M2}M_2\left(1-\frac{M_1+M_2}{K_M}\right)+\frac{\rho_2(T_S+T_R)}{a_2+T_S+T_R}-\frac{\rho_1EM_2}{a_1+E}-\mu_{M_2}M_2,\label{eq:M2}\\[2ex] \frac{dC}{dt}=&u_0-k_CC,\label{eq:C} \end{align} \end{subequations} where \begin{equation} \mu(M_2,C)=\mu_0+\mu_1M_2+\mu_2C. \label{eq:Mutation} \end{equation} The biological interpretation of the model equations is summarized below. \begin{itemize} \item The function $\mu(M_2,C)$ (\ref{eq:Mutation}) describes the rate at which drug-sensitive tumor cells acquire resistance. The parameter $\mu_0$ denotes the spontaneous mutation rate, while $\mu_1$ and $\mu_2$ quantify resistance induced by the immunosuppressive activity of M2 macrophages and chemotherapy-induced selective pressure, respectively. Furthermore, the parameter $0<\varepsilon<1$ represents the reduced sensitivity of resistant tumor cells to chemotherapy. \item Equation~(\ref{eq:TS}) describes the dynamics of drug-sensitive tumor cells. In the absence of immune responses and treatment, sensitive tumor cells are assumed to proliferate logistically with intrinsic growth rate $r_S$ and carrying capacity $K_T$, reflecting competition for nutrients and space \cite{Eftimie2021}. M2 macrophages promote tumor growth at rate $\alpha$, whereas M1 macrophages and cytotoxic CD8$^{+}$ T cells eliminate tumor cells at rates $\delta_1$ and $\kappa_1$, respectively. Chemotherapy induces tumor cell death at rate $\psi C$, while a fraction of sensitive cells acquire drug resistance according to the function $\mu(M_2,C)$. Resistant cells may revert to the drug-sensitive phenotype at rate $\eta$, representing phenotypic plasticity or re-sensitization following the removal of treatment pressure. Logistic tumor growth has been widely employed in mathematical oncology to describe growth under resource limitation \cite{Kuznetsov1994,Enderling2014}. \item Equation~(\ref{eq:TR}) describes the dynamics of drug-resistant tumor cells. Resistant tumor cells are assumed to proliferate logistically with intrinsic growth rate $r_R$ and carrying capacity $K_T$, reflecting competition for limited resources within the tumor microenvironment. Their population increases through the acquisition of resistance by drug-sensitive tumor cells at rate $\mu(M_2,C)$ and decreases through immune-mediated killing by M1 macrophages and cytotoxic CD8$^{+}$ T cells. Chemotherapy retains partial efficacy against resistant tumor cells through the term $\varepsilon\psi CT_R$, where $0<\varepsilon<1$ denotes the reduced sensitivity of resistant cells to treatment. The parameter $\eta$ represents the rate of phenotypic re-sensitization, whereby resistant tumor cells revert to the drug-sensitive phenotype following the relaxation of selective pressure. Similar sensitive--resistant tumor decompositions have been widely employed in mathematical models investigating chemotherapy resistance and tumor evolution. \item Equation~(\ref{eq:E}) describes the dynamics of cytotoxic CD8$^{+}$ T cells. Effector T cells are recruited into the tumor microenvironment at a constant rate $\Lambda_E$ and undergo tumor antigen-driven activation and clonal expansion through the saturating term \[ \frac{\rho_0(T_S+T_R)}{a_0+T_S+T_R}E, \] where $\rho_0$ denotes the maximum activation (or proliferation) rate and $a_0$ is the corresponding half-saturation constant. This functional response reflects the limited availability of tumor-associated antigens and the finite capacity of the immune system to activate effector T cells as the tumor burden increases. Cytotoxic T cells undergo natural death at rate $\mu_E$ and are further suppressed by M2 macrophages through the term $\phi M_2E$, which represents the collective immunosuppressive effects of tumor-associated macrophages, including the secretion of anti-inflammatory cytokines (e.g., IL-10 and TGF-$\beta$) and the promotion of immune checkpoint signaling that impair CD8$^{+}$ T-cell activation, proliferation and cytotoxic function \cite{Mantovani2008,Noy2014}. The parameter $\theta$ denotes density-dependent self-regulation (activation-induced cell death) of CD8$^{+}$ T-cell. \item Equation~(\ref{eq:M1}) describes the dynamics of macrophages with a dominant M1 phenotype. M1 macrophages undergo density-dependent local proliferation, modeled by a logistic growth term with carrying capacity $K_M$, representing the limited resources available for macrophage expansion within the tumor microenvironment. Activated CD8$^{+}$ T cells promote the repolarization of M2 macrophages into the M1 phenotype through the saturating term \[ \frac{\rho_1EM_2}{a_1+E}, \] where $\rho_1$ denotes the maximum repolarization rate and $a_1$ is the corresponding half-saturation constant. This term reflects the ability of activated CD8$^{+}$ T cells to secrete pro-inflammatory cytokines, particularly interferon-$\gamma$, which favor classical macrophage activation. M1 macrophages undergo natural death at rate $\mu_{M_1}$ and are depleted by tumor-mediated immunosuppressive mechanisms through the term $\delta(T_S+T_R)M_1$, where the parameter $\delta$ quantifies the inhibitory effect of increasing tumor burden on the survival and maintenance of the M1 macrophage population \cite{Mantovani2008,Murray2011}. \item Equation~(\ref{eq:M2}) describes the dynamics of macrophages with a dominant M2 phenotype. M2 macrophages undergo density-dependent local proliferation, modeled by a logistic growth term with carrying capacity $K_M$, representing the limited resources available for macrophage expansion within the tumor microenvironment. In addition, tumor cells promote the recruitment and polarization of macrophages toward the M2 phenotype through the saturating term \[ \frac{\rho_2(T_S+T_R)}{a_2+T_S+T_R}, \] where $\rho_2$ denotes the maximum tumor-induced recruitment and polarization rate, and $a_2$ is the corresponding half-saturation constant. This functional response reflects the finite capacity of tumor-derived chemokines and cytokines to recruit monocytes and promote alternative macrophage activation as tumor burden increases. The M2 macrophage population decreases through CD8$^{+}$ T-cell-mediated repolarization into the M1 phenotype and through natural mortality at rate $\mu_{M_2}$. M2 macrophages are widely recognized as key regulators of tumor progression by promoting angiogenesis, suppressing anti-tumor immune responses, facilitating extracellular matrix remodeling, and enhancing tumor invasion and metastasis \cite{Mantovani2008,Noy2014}. \item Equation~(\ref{eq:C}) describes the pharmacokinetics of chemotherapy. The drug is administered continuously at a constant infusion rate $u_0$ and is cleared from the body according to first-order kinetics with elimination rate $k_C$. This one-compartment pharmacokinetic model provides a simple autonomous representation of chemotherapy dynamics suitable for equilibrium and stability analyses. \end{itemize} Since the chemotherapy concentration satisfies \[\frac{dC}{dt}=u_0-k_CC,\] its unique equilibrium is given by \[ C_{ss}=\frac{u_0}{k_C}. \] The rapid transient dynamics of the drug concentration $C(t)$ decay almost instantaneously relative to the slow cellular state variables. Consequently, $C(t)$ rapidly equilibrates to its quasi-steady-state value, $C^*$, allowing us to decouple the fast pharmacokinetic equation from the system and evaluate the drug concentration as a constant parameter $C(t) \approx C^*$ in the slower tumor–immune dynamic equations. Consequently, the six-dimensional model can be reduced to the following five-dimensional autonomous system: \begin{subequations} \label{eq:ReducedModel} \begin{align} \frac{dT_S}{dt} =&\,r_ST_S\left(1-\frac{T_S+T_R}{K_T}\right)(1+\alpha M_2)-(\delta_1M_1+\kappa_1E)T_S-\psi C^{*}T_S-\mu(M_2,C^{*})T_S +\eta T_R,\label{eq:RTS}\\[2ex] \frac{dT_R}{dt}=&\,r_RT_R\left(1-\frac{T_S+T_R}{K_T}\right)(1+\alpha M_2)-(\delta_2M_1+\kappa_2E)T_R +\mu(M_2,C^{*})T_S-(\eta+\varepsilon\psi C^{*})T_R,\label{eq:RTR} \\[2ex] \frac{dE}{dt}=&\,\Lambda_E+\frac{\rho_0(T_S+T_R)}{a_0+T_S+T_R}E-\mu_EE-\phi M_2E-\theta E^2,\label{eq:RE}\\[2ex] \frac{dM_1}{dt} =&\,p_{M1}M_1\left(1-\frac{M_1+M_2}{K_M}\right) +\frac{\rho_1EM_2}{a_1+E} -\mu_{M_1}M_1 -\delta(T_S+T_R)M_1, \label{eq:RM1} \\[2ex] \frac{dM_2}{dt} =&\,p_{M2}M_2\left(1-\frac{M_1+M_2}{K_M}\right) +\frac{\rho_2(T_S+T_R)}{a_2+T_S+T_R} -\frac{\rho_1EM_2}{a_1+E} -\mu_{M_2}M_2, \label{eq:RM2} \end{align} \label{model} \end{subequations} where \[ C_{ss}=\frac{u_0}{k_C}, \] and \[ \mu(M_2,C^{*}) = \mu_0+\mu_1M_2+\mu_2C^{*}. \] Model parameters, their biological interpretations, baseline values, and plausible ranges to be used in the numerical simulations are presented in Table.(\ref{tab:parameters}). The parameter values were selected from published mathematical oncology and lung cancer modelling studies whenever available. Parameters for which no experimental or clinical estimates exist were assigned biologically plausible values consistent with related tumour--immune interaction models. \begin{table}[htbp] \centering \caption{Model parameters, biological interpretation, baseline values and ranges used in the simulations.} \label{tab:parameters} \footnotesize \renewcommand{\arraystretch}{1.2} \begin{tabular}{lllll} \hline Parameter & Biological meaning & Baseline & Range & References\\ \hline $r_S$ &Growth rate of drug-sensitive tumour cells &0.30 day$^{-1}$ &0.10--0.60 &\cite{Kuznetsov1994,Eftimie2021,Enderling2014}\\ $r_R$ &Growth rate of resistant tumour cells &0.24 day$^{-1}$ &0.08--0.50 &\cite{Foo2012,Gatenby2009}\\ $K_T$ &Tumour carrying capacity &$10^{9}$ cells &$10^{8}$--$10^{10}$ &\cite{Enderling2014}\\ $\alpha$ &M2-induced tumour promotion rate &$2\times10^{-8}$ day$^{-1}$cell$^{-1}$ &$10^{-9}$--$10^{-7}$ &\cite{Mantovani2008,Noy2014}\\ $\delta_1$ &M1-mediated killing of sensitive cells &$5\times10^{-8}$ day$^{-1}$cell$^{-1}$ &$10^{-9}$--$10^{-7}$ &\cite{Eftimie2021}\\ $\delta_2$ &M1-mediated killing of resistant cells &$4\times10^{-8}$ day$^{-1}$cell$^{-1}$ &$10^{-9}$--$10^{-7}$ &\cite{Eftimie2021}\\ $\kappa_1$ &CD8$^{+}$ T-cell killing of sensitive cells &$10^{-7}$ day$^{-1}$cell$^{-1}$ &$10^{-8}$--$10^{-6}$ &\cite{Kuznetsov1994,Eftimie2021}\\ $\kappa_2$ &CD8$^{+}$ T-cell killing of resistant cells &$8\times10^{-8}$ day$^{-1}$cell$^{-1}$ &$10^{-8}$--$10^{-6}$ &\cite{Kuznetsov1994,Eftimie2021}\\ $\psi$ &Chemotherapy cytotoxicity coefficient &0.40 day$^{-1}$ &0.10--1.00 &\cite{Foo2012}\\ $\varepsilon$ &Relative drug sensitivity of resistant cells &0.20 &0.05--0.50 &\cite{Foo2012,Gatenby2009}\\ $\mu_0$ &Spontaneous resistance mutation rate &$10^{-6}$ day$^{-1}$ &$10^{-8}$--$10^{-4}$ &\cite{Foo2012,Coldman1983}\\ $\mu_1$ &M2-mediated resistance induction &$10^{-8}$ day$^{-1}$cell$^{-1}$ &$10^{-9}$--$10^{-7}$ &Estimated\\ $\mu_2$ &Chemotherapy-induced resistance coefficient &0.01 day$^{-1}$ &0.001--0.05 &\cite{Foo2012}\\ $\eta$ &Re-sensitization rate &0.005 day$^{-1}$ &0--0.05 &Estimated\\ $\Lambda_E$ &Recruitment rate of CD8$^{+}$ T cells &$10^{5}$ cells day$^{-1}$ &$10^{4}$--$10^{6}$ &\cite{Kuznetsov1994,Eftimie2021}\\ $\rho_0$ &Maximum T-cell proliferation rate &0.60 day$^{-1}$ &0.20--1.00 &\cite{Eftimie2021}\\ $a_0$ &Half-saturation constant for T-cell activation &$10^{6}$ cells &$10^{5}$--$10^{8}$ &\cite{Eftimie2021}\\ $\mu_E$ &Natural death rate of CD8$^{+}$ T cells &0.05 day$^{-1}$ &0.01--0.10 &\cite{Kuznetsov1994}\\ $\phi$ & $M2$-mediated immune suppression &$10^{-8}$ day$^{-1}$cell$^{-1}$ &$10^{-9}$--$10^{-7}$ &\cite{Mantovani2008,Noy2014}\\ $p_{M1}$ &M1 proliferation rate &0.08 day$^{-1}$ &0.02--0.20 &Estimated\\ $p_{M2}$ &M2 proliferation rate &0.10 day$^{-1}$ &0.02--0.25 &Estimated\\ $K_M$ &Macrophage carrying capacity &$10^{8}$ cells &$10^{7}$--$10^{9}$ &Estimated\\ $\rho_1$ &Maximum M2$\rightarrow$M1 repolarization rate &0.20 day$^{-1}$ &0.05--0.50 &\cite{Mantovani2008,Murray2011}\\ $a_1$ &Half-saturation constant for repolarization &$10^{5}$ cells &$10^{4}$--$10^{7}$ &Estimated\\ $\mu_{M_1}$ &Natural death rate of M1 macrophages &0.03 day$^{-1}$ &0.01--0.10 &\cite{Murray2011}\\ $\delta$ &Tumour-mediated depletion of M1 macrophages &$10^{-8}$ day$^{-1}$cell$^{-1}$ &$10^{-9}$--$10^{-7}$ & Estimated\\ $\rho_2$ &Maximum recruitment of M2 macrophages &$5\times10^{5}$ cells day$^{-1}$ &$10^{4}$--$10^{6}$ &\cite{Mantovani2008,Noy2014}\\ $a_2$ &Half-saturation constant for M2 recruitment &$10^{6}$ cells &$10^{5}$--$10^{8}$ &Estimated\\ $\mu_{M_2}$ &Natural death rate of $M2$ macrophages &0.03 day$^{-1}$ &0.01--0.10 &\cite{Murray2011}\\ $u_0$ &Drug infusion rate &1.0 concentration units day$^{-1}$ &0.2--5.0 &\cite{dePillis2005}\\ $k_C$ &Drug elimination rate &0.50 day$^{-1}$ &0.20--1.50 &\cite{dePillis2005}\\ $\theta$ &Density-dependent self-regulation of CD8+ T cells &$5\cdot 10^{-7}$ day$^{-1}$ &$10^{-8}$--$10^{-6}$&Estimated\\ \hline \end{tabular} \end{table} \subsection{Positivity and boundedness of solutions} Before investigating the qualitative dynamics of system \eqref{eq:ReducedModel}, we establish that it is mathematically and biologically well posed by showing that solutions remain non-negative and uniformly bounded for all future time. \begin{theorem} Let \[ (T_S(0),T_R(0),E(0),M_1(0),M_2(0)) \in \mathbb{R}_+^5. \] Then the unique solution of system \eqref{eq:ReducedModel} exists for all \(t\ge0\), remains non-negative, and is uniformly bounded in a compact positively invariant set of \(\mathbb{R}_+^5\). \end{theorem} \begin{proof} Since the right-hand side of system \eqref{eq:ReducedModel} is continuously differentiable, existence and uniqueness of solutions follow from the Picard--Lindelöf theorem \cite{Coddington1955,Perko2001}. To establish positivity, suppose one of the state variables reaches zero for the first time. Then \begin{align*} \left.\frac{dT_S}{dt}\right|_{T_S=0}&=\eta T_R\ge0,\qquad \left.\frac{dT_R}{dt}\right|_{T_R=0}=\mu(M_2,C^*)T_S\ge0,\qquad \left.\frac{dE}{dt}\right|_{E=0}=\Lambda_E>0,\cr \left.\frac{dM_1}{dt}\right|_{M_1=0}&=\frac{\rho_1EM_2}{a_1+E}\ge0,\qquad \left.\frac{dM_2}{dt}\right|_{M_2=0}=\frac{\rho_2(T_S+T_R)}{a_2+T_S+T_R} \ge0. \end{align*} Hence the vector field is inward-pointing (or tangent) on each coordinate hyperplane, implying that \(\mathbb{R}_+^5\) is positively invariant. Next, define the total tumour and macrophage populations as \[ T=T_S+T_R, \qquad M=M_1+M_2. \] Adding equations \eqref{eq:RM1} and \eqref{eq:RM2} yields \[ \frac{dM}{dt} \le p_{\max}M -\frac{p_{\min}}{K_M}M^2 +\rho_2, \] where $p_{\max}=\max\{p_{M1},p_{M2}\},$ and $p_{\min}=\min\{p_{M1},p_{M2}\}.$ By the comparison theorem, $M(t)\le M^*,$ where \[ M^* = \frac{K_M}{2p_{\min}} \left( p_{\max} + \sqrt{ p_{\max}^2+ \frac{4p_{\min}\rho_2}{K_M}} \right). \] Consequently, $M_1(t),\,M_2(t)\le M^*.$ Similarly, adding equations \eqref{eq:RTS} and \eqref{eq:RTR} gives \[ \frac{dT}{dt} \le r_{\max} \left(1-\frac{T}{K_T}\right) (1+\alpha M^*)T, \] where $r_{\max}=\max\{r_S,r_R\}.$ Comparison with the logistic equation implies $T(t)\le K_T.$ Finally, from equation \eqref{eq:RE}, \[ \frac{dE}{dt} \le \Lambda_E + (\rho_0-\mu_E)E -\theta E^2. \] The corresponding comparison equation possesses the positive equilibrium \[ E^* = \frac{ (\rho_0-\mu_E) + \sqrt{ (\rho_0-\mu_E)^2 + 4\theta\Lambda_E}} {2\theta}, \] and therefore $E(t)\le E^*.$ Hence every solution eventually enters the compact region \[ \Omega= \left\{ (T_S,T_R,E,M_1,M_2)\in\mathbb{R}_+^5: \; T\le K_T,\; E\le E^*,\; M\le M^* \right\}, \] which is positively invariant. Therefore, all solutions of system \eqref{eq:ReducedModel} remain non-negative and uniformly bounded for all \(t\ge0\). \end{proof} \subsection{Model steady states and their stability} In this section, we conduct a thorough analysis of the steady states of the model and evaluate their stability. Understanding these steady states is crucial, as they provide insights into the long-term behavior of the system. By examining the conditions under which these steady states are maintained or disrupted, we gain valuable information that can inform the design of effective interventions. This detailed exploration allows us to identify the factors that influence stability and to develop strategies that can promote desirable outcomes within the system \cite{Noy2014}. This analysis begins with the mathematical model articulated in Eq.(\ref{eq:TR}) through Eq(\ref{eq:C}), which represents a comprehensive six-dimensional ordinary differential equation (ODE) system.\\ This system describes the complex interactions and dynamics of various biological components involved in tumor progression and treatment response. Specifically, it includes sensitive tumor cells ($T_{S}$), which are vulnerable to therapeutic interventions; resistant tumor cells ($T$), capable of evading treatment; effector immune cells ($E$), which play a crucial role in targeting and eliminating tumor cells; M1 macrophages ($M_{1}$), known for their pro-inflammatory and tumoricidal properties; M2 macrophages ($M_{2}$), which are associated with tumor promotion and immune suppression; and the concentration of the drug ($C$) being administered \cite{Eftimie2021}. The local stability analysis is conducted by examining the linear stability around the steady states of this multidimensional system. By doing so, we can systematically assess how small perturbations in the system may affect its overall dynamics. This analysis not only enhances our understanding of the intrinsic behavior of the tumor-immune dynamics but also aids in predicting potential outcomes of different treatment strategies, ultimately contributing to the development of more effective therapeutic approaches \cite{Mantovani2008}.\\ Let us define the vector $X$ as $ X = (T_{S}, T_{R}, E, M_{1}, M_{2}, C)^{T} $, where each component represents a specific variable in our system. From this definition, we can construct a $6 \times 6 $ Jacobian matrix, denoted $ J(X) $, which is composed of partial derivatives of each variable in $X$ with respect to another variable. Specifically, the matrix $J(X)$ is given by $ J(X) = \left[\frac{\partial x_{i}}{\partial y_{j}}\right] $. This Jacobian encapsulates the relationships between the variables defined in Eq.(\ref{eq:TS}) through in Eq.(\ref{eq:C}), resulting in the following matrix: \[ J= \left( \begin{tabular}{llllll} $J_{11}$ & $J_{12}$ & $J_{13}$ & $J_{14}$ & $J_{15}$ & $J_{16}$\\ $J_{21}$ & $J_{22}$ & $J_{23}$ & $J_{24}$ & $J_{25}$ & $J_{26}$\\ $J_{31}$ & $J_{32}$ & $J_{33}$ & $0$ & $J_{35}$ & $0$\\ $J_{41}$ & $J_{42}$ & $J_{43}$ & $J_{44}$ & $J_{45}$ & $0$\\ $J_{51}$ & $J_{52}$ & $J_{53}$ & $J_{54}$ & $J_{55}$ & $0$\\ $0$ & $0$ & $0$ & $0$ & $0$ & $-k_{C}$\\ \end{tabular} \right) \] where the non-zero elements are grouped into two categories which are tumor subsystem rows and immune and microenvironment rows as follows \begin{itemize} \item[(i)]Tumor Subsystem Rows ( $T_{S},T_{R}$) \begin{eqnarray} \ds J_{11}&=&r_{S}\left(1-\frac{2T_{S}+T_{R}}{K_{T}}\right)+\alpha M_{2}-\delta _{1}M_{1}-K_{1}E-\psi C-\mu (M_{2},C)\label{eq1}\\[5pt] \ds J_{12}&=&\frac{-r_{S}T_{S}}{K_{T}}+\eta\label{eq2}\\[5pt] \ds J_{13}&=&-k_{1}T_{S},~~J_{14}=-\delta _{1}T_{S},~~J_{15}=(\alpha -\mu _{1})T_{S},~~J_{16}=-(\psi +\mu _{2})T_{S}\label{eq3}\\[5pt] \ds J_{21}&=&-\frac{r_{R}T_{R}}{K_{T}}+\mu (M_{2}, C)\label{eq4}\\[5pt] \ds J_{22}&=&r_{R}\left(1-\frac{T_{s}+2T_{R}}{K_{T}}\right)+\alpha M_{2}-\delta_{2} M_{1}-k_{2}E-\eta -\varepsilon \psi C\label{eq5}\\[5pt] \ds J_{23}&=&-k_{2}T_{R},~~J_{24}=-\delta_{2} T_{R},~~J_{25}=\alpha T_{R}+\mu _{1}T_{S},~~J_{26}=\mu _{2}T_{S}-\varepsilon \psi T_{R}~~~~~~\label{eq6} \end{eqnarray} \item[(ii)] Immune and Microenvironment Rows ($E,M_{1},M_{2}$) \begin{eqnarray} \ds J_{31}&=&J_{32}=\frac{\rho _{0}a_{0}E}{\left(a_{0}+T_{S}+T_{R}\right)^{2}}\label{eq7}\\[5pt] \ds J_{33}&=&\frac{\rho _{0}(T_{S}+T_{R})}{a_{0}+T_{S}+T_{R}}-\mu E-\phi M_{2}-\theta E,~~J_{35}=-\phi E\label{eq8}\\[5pt] \ds J_{41}&=&J_{42}=-\delta M_{1},~~J_{43}=\frac{\rho _{1}a_{1}M_{2}}{\left(a_{1}+E\right)^{2}}\label{eq9}\\[5pt] \ds J_{44}&=&p_{M_{1}}\left(1-\frac{2M_{1}+M_{2}}{K_{M}}\right)-\mu _{M_{1}}-\delta(T_{S}+T_{R})\label{eq10}\\[5pt] \ds J_{45}&=&-\frac{p_{M_{1}}M_{1}}{K_{M}}+\frac{\rho _{1}E}{a_{1}+E}\label{eq11}\\[5pt] \ds J_{51}&=&J_{52}=\frac{\rho _{2}a_{2}}{\left(a_{2}+T_{S}+T_{R}\right)^{2}},~~J_{53}=-\frac{\rho _{1}a_{1}M_{2}}{\left(a_{1}+E\right)^{2}}\label{eq12}\\[5pt] \ds J_{54}&=&-\frac{p_{M_{2}}M_{2}}{K_{M}},~~J_{55}=p_{M_{2}}\left(1-\frac{M_{1}+2M_{2}}{K_{M}}\right)-\frac{\rho _{1}E}{a_{1}+E}-\mu _{M_{2}}~~~~~~~~~~~~~~~\label{eq13} \end{eqnarray} \end{itemize} In the complete absence of tumor populations ($T_{S} = 0, T_{R} = 0$) and macrophage lineages ($M_{1} = 0, M_{2} = 0$), along with zero remaining drug concentration ($C = 0$), the system reduces to its uninfected and uninflamed state. Under these conditions, the system admits a disease-free (or tumor-free) steady state, denoted by $\hat{E}^{*}$, where the effector cell population $E$ equilibrates at its constant physiological baseline level determined strictly by its baseline recruitment and natural death rate. The tumor steady state of Eq.(\ref{eq:TR}) through Eq.(\ref{eq:C}) is given by \[\hat{E}^{*}=(0,0,E_{0},0,0,C^{*})\] where $C^{*}=\frac{u_{0}}{k_{C}}$ and $E_{0}$ is the unique positive solution to $\Lambda _{E}-\mu _{E}E-\theta E^{2}=0$, then \[E_{0} = \frac{-\mu_E + \sqrt{\mu_E^2 + 4\theta \Lambda_E}}{2\theta}\] At $\varepsilon _{0}$, off-diagonal blocks decouple the system into smaller sub-blocks. The Jacobian simplifies to: \[ J(\varepsilon _{0})= \left( \begin{tabular}{llll} $J_{T}$ & $0$ & $0$ & $0$ \\ $J_{ET}$ & $\lambda _{E}$ & $J_{ME}$ & $0$\\ $0$ & $0$ & $J_{M}$ & $0$\\ $0$ & $0$ & $0$ & $-k_{C}$ \\ \end{tabular} \right) \] giving Chemotherapy and Effector eigenvalues \begin{eqnarray} \ds \lambda _{1}&=&-k_{C}< 0 \label{eq14}\\[5pt] \ds \lambda _{2}&=&\frac{\partial f_{E}}{\partial E}|_{\varepsilon _{0}}=-\mu _{E}-2\theta E_{0}=-\mu _{E}^{2}+4\theta \Lambda _{E}<0 \label{eq15} \end{eqnarray} Furthermore, for macrophage subsystem eigenvalues we introduced the submatrix $J_{M}$ governing macrophage invasion, which is upper triangular, then the eigenvalues are the diagonal elements given as follows: \begin{eqnarray} \ds \lambda _{M_{1}}&=&p_{M_{1}}-\mu _{M_{1}},~~\lambda _{M_{2}}=p_{M_{2}}-\mu _{M_{2}}-\frac{\rho _{1}E_{0}}{a_{1}+E_{0}}\label{16} \end{eqnarray} To guaranty the local asymptotic stability of the tumor-free equilibrium state $\mathcal{E}_0$ against macrophage infiltration, the parameters of the system must satisfy two specific analytical threshold conditions sating that \(p_{M_1} < \mu_{M_1} ~~\text{and}~~ p_{M_2} < \mu_{M_2} - \frac{\rho_1 E_0}{a_1 + E_0}\). Analogously to the reduced system presented in Eq.(\ref{eq:RTS}) to Eq.(\ref{eq:RM2}), operating under the assumption that chemotherapy pharmacokinetics equilibrate on a significantly faster timescale than the underlying tumor–immune interactions, maintaining mathematical stability ensures that the dynamical system yields consistent, structurally predictable behaviors across varying initial conditions and parameter inputs \cite{miller2018stable}. \subsubsection{Reproduction number of tumor cells population} The reproductive number for the system is crucial in the stability and steady analysis derived from the next-generation matrix technique that represent steady-states the spread rate of cancer depending on the value of $R_{0}$. After substituting the value of parameters we get that the system is disease free due to introducing new control variables. Then \begin{eqnarray} R_0 = \frac{(\beta_S d_R + \beta_R d_S) + \sqrt{(\beta_S d_R - \beta_R d_S)^2 + 4 \beta_S \beta_R \eta \mu^*}}{2(d_S d_R - \eta \mu^*)} \end{eqnarray} where the total effective clearance or loss rates $d_S$ and $d_R$ at $\hat{E}^*$ are defined as $d_S = \delta_{1} M_{1}^{*} + \kappa_1 E^{*} + \psi C^{*} + \mu^{*}$ and $d_{R} = \delta_{2} M_{1}^{*} + \kappa_{2} E^{*} + \varepsilon \psi C^{*} + \eta$, see Appendix.(\ref{app:extra2}) for full working. \begin{theorem} Assuming that $R_0 > 1$ such that a unique positive endemic equilibrium given by: \[\hat{E}^* = (T_S^*, T_R^*, E^*, M_1^*, M_2^*, C^*)\] exists in the interior of the positive invariant region $\Omega \subset \mathbb{R}_+^6$. Then $\hat{E}^*$ is Globally Asymptotically Stable (GAS) in $\mathrm{Int}(\Omega)$ if the inter-compartmental interaction terms satisfy the parameter constraints bounded by matrix $Q \mathbf{z} \le 0$ where $Q$ is real symmetric and negative semi-definite. \label{theorem1.3} \end{theorem} \begin{proof} Defining the equilibrium relations at $\hat{E}^*$, all time derivatives in Eq.(\ref{eq:TR}) through Eq.(\ref{eq:C}) are zero ($\frac{d T_S}{dt} = \frac{d T_R}{dt} = \frac{d E}{dt} = \frac{d M_1}{dt} = \frac{d M_2}{dt} = \frac{d C}{dt} = 0$). Dividing each differential equation by its respective state variable evaluated at $E^*$ yields the baseline identity balance equations: \[\begin{aligned} 0 &= r_S\left(1-\frac{T_S^*+T_R^*}{K_T}\right) + \alpha M_2^* - \delta_1 M_1^* - \kappa_1 E^* - \psi C^* - \mu^* + \eta \frac{T_R^*}{T_S^*} \\ 0 &= r_R\left(1-\frac{T_S^*+T_R^*}{K_T}\right) + \alpha M_2^* - \delta_2 M_1^* - \kappa_2 E^* + \mu^*\frac{T_S^*}{T_R^*} - \eta - \varepsilon \psi C^* \\ 0 &= \frac{\Lambda_E}{E^*} + \frac{\rho_0(T_S^*+T_R^*)}{a_0+T_S^*+T_R^*} - \mu_E - \phi M_2^* - \theta E^* \\ 0 &= p_{M1}\left(1-\frac{M_1^*+M_2^*}{K_M}\right) + \frac{\rho_1 E^* M_2^*}{M_1^* (a_1+E^*)} - \mu_{M1} - \delta(T_S^*+T_R^*) \\ 0 &= p_{M2}\left(1-\frac{M_1^*+M_2^*}{K_M}\right) + \frac{\rho_2(T_S^*+T_R^*)}{M_2^* (a_2+T_S^*+T_R^*)} - \frac{\rho_1 E^*}{a_1+E^*} - \mu_{M2} \\ 0 &= \frac{u_0}{C^*} - k_C \end{aligned}\] where $\mu^* = \mu(M_2^*, C^*) = \mu_0 + \mu_1 M_2^* + \mu_2 C^*$. Following the Lyapunov stability framework established by Ahmad et al. \cite{ahmad2024mathematical}, we define the positive definite function $V: \mathrm{Int}(\Omega) \to \mathbb{R}_{\ge 0}$: \[V(T_S, T_R, E, M_1, M_2, C) = w_1 g\left(\frac{T_S}{T_S^*}\right) + w_2 g\left(\frac{T_R}{T_R^*}\right) + w_3 g\left(\frac{E}{E^*}\right) + w_4 g\left(\frac{M_1}{M_1^*}\right) + w_5 g\left(\frac{M_2}{M_2^*}\right) + w_6 g\left(\frac{C}{C^*}\right)\] where $g(z) = z - 1 - \ln z \ge 0$ for all $z > 0$, with $g(z) = 0$ if and only if $z = 1$.Taking the derivative of $V$ with respect to time $t$: \[\frac{dV}{dt} = w_1 \left(1 - \frac{T_S^*}{T_S}\right) \frac{dT_S}{dt} + w_2 \left(1 - \frac{T_R^*}{T_R}\right) \frac{dT_R}{dt} + w_3 \left(1 - \frac{E^*}{E}\right) \frac{dE}{dt}\] \[ +w_4 \left(1 - \frac{M_1^*}{M_1}\right) \frac{dM_1}{dt} + w_5 \left(1 - \frac{M_2^*}{M_2}\right) \frac{dM_2}{dt} + w_6 \left(1 - \frac{C^*}{C}\right) \frac{dC}{dt}\] Decomposing $\frac{dV}{dt}$ into component terms, substituting the structural expressions of the system equations and using the equilibrium conditions, we group $\frac{dV}{dt}$ into two distinct parts which are the quadratic density-dependent terms ($\mathcal{D}$) and non-linear interaction or coupling terms ($\mathcal{I}$): \[\frac{dV}{dt} = \mathcal{D} + \mathcal{I}\] 1. Quadratic Density-Dependent terms ($\mathcal{D}$) representing intrinsic logistic self-limitation, natural mortality, decay, and crowding: \[\begin{aligned} \mathcal{D} = &- w_1 \frac{r_S}{K_T} (T_S - T_S^*)^2 - w_2 \frac{r_R}{K_T} (T_R - T_R^*)^2 - w_3 \theta (E - E^*)^2 \\ &- w_4 \frac{p_{M1}}{K_M} (M_1 - M_1^*)^2 - w_5 \frac{p_{M2}}{K_M} (M_2 - M_2^*)^2 - w_6 k_C \frac{(C - C^*)^2}{C} \end{aligned}\] Since $w_i, K_T, K_M, \theta, k_C, r_S, r_R > 0$, every term in $\mathcal{D}$ is strictly negative for all $\mathbf{x} \neq \mathbf{x}^*$. \\ 2. Cross coupling and saturation terms ($\mathcal{I}$) using the standard identity $1 - z + \ln z \le 0$ for Volterra expressions, the non-linear saturation terms are divided into bounded Volterra structures and cross-deviations: \\ Saturating Immune Recruitment Terms: \[w_3 \left(1 - \frac{E^*}{E}\right) \left[ \frac{\rho_0 (T_S + T_R)}{a_0 + T_S + T_R} E - \frac{\rho_0 (T_S^* + T_R^*)}{a_0 + T_S^* + T_R^*} E \right] \le w_3 \frac{\rho_0 a_0 (E - E^*)(T_S - T_S^* + T_R - T_R^*)}{(a_0 + T_S + T_R)(a_0 + T_S^* + T_R^*)}\] Inter-tumor Switching or Mutation and Resensitization: \[w_1 \eta \left( 1 - \frac{T_S^*}{T_S} \right) \left( T_R - T_R^* \frac{T_S}{T_S^*} \right) + w_2 \mu^* \left( 1 - \frac{T_R^*}{T_R} \right) \left( T_S - T_S^* \frac{T_R}{T_R^*} \right)\] Using the inequality $(2 - \frac{x}{x^*} - \frac{x^*}{x}) \le 0$, these terms evaluate to non-positive values when balanced by setting $w_1 \eta T_R^* = w_2 \mu^* T_S^*$. Quadratic Form Bounds and Sufficient Conditions. Let \[\mathbf{z} = \begin{pmatrix} \vert{}T_S - T_S^*\vert{}, & \vert{}T_R - T_R^*\vert{}, & \vert{}E - E^*\vert{}, & \vert{}M_1 - M_1^*\vert{}, & \vert{}M_2 - M_2^*\vert{}, & \vert{}C - C^*\vert{} \end{pmatrix}^T\] By applying Young's inequality ($2xy \le \epsilon x^2 + \frac{1}{\epsilon} y^2$) to bound the off-diagonal coupling terms in $\mathcal{I}$ using the negative diagonal terms in $\mathcal{D}$, we can express $\frac{dV}{dt}$ as a matrix quadratic form: \[\frac{dV}{dt} \le -\mathbf{z}^T Q \mathbf{z}\] where $Q = (Q_{ij}) \in \mathbb{R}^{6 \times 6}$ is a symmetric matrix ($Q = Q^T$) with diagonal elements $Q_{ii}$ for $i = 1, \dots, 6$, thus \[\begin{aligned} Q_{11} &= w_1 \frac{r_S}{K_T} - \epsilon_1, \quad Q_{22} = w_2 \frac{r_R}{K_T} - \epsilon_2, \quad Q_{33} = w_3 \theta - \epsilon_3 \\ Q_{44} &= w_4 \frac{p_{M1}}{K_M} - \epsilon_4, \quad Q_{55} = w_5 \frac{p_{M2}}{K_M} - \epsilon_5, \quad Q_{66} = w_6 k_C - \epsilon_6 \end{aligned}\] Where, $\epsilon_i > 0$ depend on the upper bounds of the interaction parameters ($\delta_1, \kappa_1, \alpha, \phi, \psi$). Provided the intra-specific competition or death rates exceed the inter-specific interaction rates ($Q$ is positive definite, meaning all principal minors of $Q$ are strictly positive), we have: \[\frac{dV}{dt} \le 0 \quad \forall \mathbf{x} \in \mathrm{Int}(\Omega)\] Applying the LaSalle's invariance principle negative semi-definiteness: $\frac{dV}{dt} \le 0$ for all $(T_S, T_R, E, M_1, M_2, C) \in \mathrm{Int}(\Omega)$.Set of Zero Derivatives: Let $M = \left\{ \mathbf{x} \in \mathrm{Int}(\Omega) \;\Big\vert{}\; \frac{dV}{dt} = 0 \right\}$.$\frac{dV}{dt} = 0$ holds if and only if each quadratic difference vanishes: \[(T_S - T_S^*)^2 = 0, \quad (T_R - T_R^*)^2 = 0, \quad (E - E^*)^2 = 0, \quad (M_1 - M_1^*)^2 = 0, \quad (M_2 - M_2^*)^2 = 0, \quad (C - C^*)^2 = 0\] This implies $M = \{ \hat{E}^* \}$. By LaSalle's Invariance Principle, every solution trajectory starting in $\mathrm{Int}(\Omega)$ approaches $\hat{E}^*$ as $t \to \infty$.Thus, the endemic equilibrium $\hat{E}^*$ is Globally Asymptotically Stable (GAS). \end{proof} \section{Numerical Simulations} In this section, we establish the baseline parameter values and initial state vector $x(0)$ given by \[\mathbf{x}(0) = [T_S(0), T_R(0), E(0), M_1(0), M_2(0), C(0)]^T\] To rigorously characterize the competitive dynamics governing lung cancer progression, this study formulates a deterministic mathematical framework aimed at conducting quantitative parameter sensitivity analysis and identifying the dominant kinetic drivers of tumor growth, treatment response, and immune-mediated clearance \cite{elinav2013inflammation}. To capture this early phase of anti-tumor immunosurveillance prior to the onset of profound tumor-induced immunosuppression, we specify positive, non-zero baseline recruitment and intrinsic growth rates for both cytotoxic effector lymphocytes ($E$) and tumoricidal $M_1$-polarized macrophages ($M_1$). Modeling these baseline populations explicitly ensures that the host immune microenvironment actively engages the nascent tumor load at inception, establishing a physiological baseline against which subsequent immune exhaustion, phenotypic switching, and therapeutic interventions can be systematically evaluated. To thoroughly evaluate the structural stability and robustness of the formulated ordinary differential equation (ODE) framework, we execute a bifurcated sensitivity analysis encompassing both local and global methodologies. Local Sensitivity Analysis (LSA) focuses on computing normalized sensitivity indices (elasticities) via direct partial derivatives evaluated at the baseline parameter set. This quantifies how localized, small-scale perturbations in individual parameter values ($p_i$) immediately impact state outputs and thresholds such as the basic reproduction number ($R_0$) \cite{donze2011robustness}. Global Sensitivity Analysis (GSA) utilizes multi-dimensional sampling techniques—such as Latin Hypercube Sampling paired with Partial Rank Correlation Coefficients (LHS-PRCC) or variance based Sobol' indices across the entire feasible parameter space \cite{donze2011robustness, elinav2013inflammation}. \subsection{Baseline System Dynamics} Baseline models capture the natural progression of cancer and immune system interactions, tumor growth kinetics, and cellular population dynamics over time. As illustrated in Fig.(\ref{Fig:1 Sensitive model simulation}), the numerical simulation captures the temporal dynamics of each compartment within the tumor–immune chemotherapy system. Beginning with the sensitive tumor cell population ($T_S$), the profile exhibits an immediate and rapid decline from an initial burden of $10^6$ cells. This acute phase reflects the primary cytotoxic effect of chemotherapy ($C$), coupled with the targeted anti-tumor activities of cytotoxic macrophages M1 ($M_1$) and effector immune cells ($E$) \cite{Eftimie2021}. Mechanically, during this early phase, combined elimination rates comprising chemotherapeutic clearance effector-mediated lysis, and macrophage phagocytosis of M1 substantially outperform intrinsic logistic renewal. During the extended $300$-day simulation period, $T_S$ undergoes continuous and sustained decay, ultimately reaching near-zero or trace levels. This asymptotic trajectory toward near-eradication indicates that under the specified parameter regime, the synergistic action of chemotherapeutic administration and endogenous immune surveillance is sufficient to suppress sensitive cell proliferation and drive the population toward extinction \cite{bracci2014immune}.\\ Resistant Tumor Cells ($T_R$) have initial growth starting at $10^3$ cells and initially experiences growth. This growth occurs even as $T_S$ cells are declining, which can be attributed to the $\mu(M_2,C)T_S$ term, representing the mutation of sensitive cells into resistant ones, and the proliferation term, which allows resistant cells to proliferate. Resistant cells are less susceptible to chemotherapy ($-\epsilon \psi CT_{R}$ where $\epsilon$ is 0.2, meaning $20\%$ sensitivity) \cite{kareva2015metronomic}. After an initial increase, the $T_R$ population reaches a peak around $100-150$ days before also starting to decline, eventually reaching very low levels. The decline in $T_R$ suggests that even though they are resistant to chemotherapy, the combined effects of $M1$ macrophages and effector T-cells eventually overcome their proliferation, especially as the total tumor burden decreases \cite{basak2023tumor}. Moving to the effector cells ($E$) population which starts at $10^{4}$ cells and exhibits a sharp increase in the early phases of the simulation, peaking around $5 \times 10^{5}$ cells. This rapid expansion is driven by the recruitment rate $\Lambda_E$ and proliferation due to the presence of tumor cells . After reaching its peak, the $E$ population remains relatively high for a significant period before slowly declining \cite{tufail2025immune}. This sustained presence is crucial for controlling both $T_S$ and $T_R$. \begin{center} \begin{figure}[hbt!] \centering %\fbox{% %\tmpframe{ \includegraphics[width = 6.3in]{image0.1.png} %\includegraphics[width = 6.0in]{image0.2.png} %} %} \caption{ \textit{ The baseline plots offer a detailed comparison of six model variables, showing their behavior under three distinct sets of initial conditions. This analysis enhances our understanding of each variable's response to different scenarios.}}. \label{Fig:1 Sensitive model simulation} \end{figure} \end{center} The eventual decline likely corresponds to the reduction in the total tumor burden, which reduces the proliferative stimulus, and the effects of natural death, $M2$ suppression, and density-dependent regulation \cite{basak2023tumor}. Furthermore the $M_1$ population starts at $10^3$ cells and shows a strong, continuous increase throughout the simulation, reaching levels well into the millions. This suggests a robust $M_{1}$ response, indicating effective repolarization from $M_{2}$ macrophages and proliferation. The high levels of $M_1$ contribute significantly to the killing of both sensitive and resistant tumor cells. The growth eventually slows down as it approaches the macrophage carrying capacity ($K_{M}$), but it remains a dominant immune component \cite{basak2023tumor}. The $M_2$ population starts at $10^4$ cells and also initially increases. This can be driven by recruitment due to tumor presence. However, after reaching a peak (around $9 \times 10^5$), the $M_2$ population starts to decline, eventually reaching very low levels. This suggests that the repolarization to $M1$ macrophages becomes dominant, or their recruitment sources diminish as the tumor burden decreases \cite{van2018molecular}.\\ The drug concentration starts at $0.1$ and rapidly increases due to the continuous infusion ($u_{0}$). It quickly reaches a steady state around $1.9-2.0$ concentration units. This steady state is determined by the balance between the infusion rate $u_{0}$ and the elimination rate ($-k_{C}C$). The sustained high concentration of $C$ ensures continuous pressure on sensitive tumor cells and contributes to the mutation of $T_{S}$ to $T_{R}$ \cite{liu2025drug}. The simulation suggests that, with the chosen parameters, the treatment strategy (chemotherapy combined with the immune response) is largely effective in controlling both sensitive and resistant tumor populations over the simulated period \cite{liu2025drug, bracci2014immune}. The initial decline of sensitive cells is robust and while resistant cells initially grow, they are eventually brought under control by the persistent immune response, particularly the strong population of macrophages $M1$ and effector T cells. The chemotherapy drug reaches and maintains a steady, effective concentration throughout. This outcome implies a successful therapeutic intervention in this specific model scenario \cite{bracci2014immune}.In order to gain further insights in the baseline system, we included AUC (Area Under the Curve) shown in Table.(\ref{tab:AUC}) which helps us quantify the cumulative tumor burden over the entire simulation period. \begin{table}[htbp] \centering \caption{Combined Comparison For Total Tumor Burden AUC} \label{tab:AUC} \footnotesize \renewcommand{\arraystretch}{1.2} \begin{tabular}{lrrr} \hline & Original IC & First New IC& Second New IC \\ Metric & & & \\ \hline $T_{S}$ (AUC) & $4.27\times 10^{6}$ & $7.01\times 10^{6}$ & $1.53\times 10^{6}$ \\ $T_{R}$ (AUC) & $1.90\times 10^{6}$ & $3.11\times 10^{6}$ & $1.53\times 10^{6}$ \\ Total Tumor Burden (AUC) & $6.17\times 10^{6}$ & $1.01\times 10^{7}$ & $3.06\times 10^{6}$ \\ \hline \end{tabular} \end{table} Instead of just looking at the tumor size at a single point in time like the steady-state value, AUC provides a single metric that integrates the tumor population across all simulated days. This means it accounts for both how large the tumor gets and for how long it stays at certain levels. Moreover, it allows for a straightforward comparison between different initial conditions or treatment strategies. The results in the table can also be illustrated with a bar chart in Fig.(\ref{Fig:1 Bar plots AUC}) explaining that lower Area Under the Curve (AUC) value is often associated with a more favorable clinical outcome. \begin{center} \begin{figure}[hbt!] \centering %\fbox{% %\tmpframe{ \includegraphics[width = 5.5in]{image2.2.png} %\includegraphics[width = 3.3in]{image0.3.png} %} %} \caption{ \textit{The bar chart provides a visual comparison of the total tumor burden against area under the curve (AUC) for the three different initial conditions (IC).}}. \label{Fig:1 Bar plots AUC} \end{figure} \end{center} This relationship suggests that throughout the treatment or observation period, there was a lower overall tumor burden. Essentially, a reduced AUC indicates that the tumor size or activity was less pronounced, which may correlate with better responses to therapy and improved patient prognosis. \subsection{Local Sensitivity Analysis} To account for the uncertainty inherent in range-estimated parameters, we investigate the robustness of the two models presented in Eq.(\ref{eq:TS}) to Eq.(\ref{eq:RTR}) by examining tumor size sensitivity to small perturbations in initial conditions and systemic parameters. Initial local sensitivity assessments evaluate the isolated effects of baseline conditions and macrophage polarization or re-polarization dynamics. To overcome the limitations of local methods and account for non-linear interactions across the parameter space, we extend this to a global sensitivity analysis \cite{ma2024comprehensive}. In order to perform a local sensitivity analysis, we will need to perturb each parameter by a small percentage ($\pm 1\%$) and observe the change in a chosen output metric. For this example, we will look at the final values of $TS$ and $TR$ at the end of the simulation. Then, we calculate the relative sensitivity index for each parameter using the formulae: \begin{eqnarray} S_P^X =\frac{P}{X} \left(\frac{\partial X}{\partial P}\right) \approx \frac{(X_{P + \Delta P} - X_P) / X_P}{\Delta P / P}\label{eq.sa} \end{eqnarray} The output with the baseline parameter value is denoted as $X_P$, while $X_{P + \Delta P}$ represents the output when the parameter $P$ is altered by $\Delta P$. Analyzing the results, we find that $\vert{}S_P^X\vert{} > 1$ indicates high sensitivity, meaning that the changes in output surpass the magnitude of the parameter perturbation. Such parameters have a significant influence on system behavior, necessitating precise estimation \cite{jacobs1990risk}. Conversely, $\vert{}S_P^X\vert{} \approx 0$ reflects low sensitivity, implying that the output remains largely unaffected by changes in the parameter, thus making their calibration less critical. A positive sensitivity ($S_P^X > 0$) indicates that an increase in parameter $P$ results in an increase in output $X$, while a negative sensitivity ($S_P^X < 0$) signifies that an increase in $P$ leads to a decrease in $X$. \begin{center} \begin{figure}[hbt!] \centering %\fbox{% %\tmpframe{ \includegraphics[width = 6.5in]{image0.2.png} %\includegraphics[width = 3.3in]{image0.3.png} %} %} \caption{ \textit{The sensitivity indices for \( T_{S} \) and \( T_{R} \) help to elucidate which parameters exert the most significant influence on their final values. By analyzing these indices, we can identify and prioritize the parameters that may require more focused attention or adjustment to achieve desired outcomes.}}. \label{Fig:1 Bar plots} \end{figure} \end{center} To gain a deeper insight into the parameters that significantly influence cell population, we conducted a local sensitivity analysis and summarized our findings using bar plots, as depicted in Fig.(\ref{Fig:1 Bar plots}) derived from Eq.(\ref{eq.sa}). In these plots, the length of each bar visually indicates the magnitude of a parameter's effect on the cell population. Specifically, longer bars reflect a greater impact, showcasing the importance of these parameters in determining cell counts. Additionally, the orientation of the bars adds further context: bars extending to the right signify that an increase in the parameter's value yields a corresponding increase in cell population (indicating positive sensitivity), while bars extending to the left indicate that an increase in the parameter leads to a decrease in cell population (denoting negative sensitivity). This visual representation makes it easier to identify parameters that merit further investigation and may be crucial in modeling the system. The parameters are typically arranged from top to bottom (or vice versa) based on the absolute magnitude of their sensitivity, enabling a quick identification of the most critical factors. \begin{center} \begin{figure}[hbt!] \centering %\fbox{% %\tmpframe{ \includegraphics[width = 6.3in]{image0.62.png} %\includegraphics[width = 3.3in]{image0.61.png} %} %} \caption{ \textit{The local sensitivity analysis for $K_M$ and $\kappa_1$ on the final $T_S$ and $T_R$ populations.}}. \label{Fig:1 actual tumor growth} \end{figure} \end{center} To provide a deeper understanding of the local sensitivity analysis, we present visual representations in Fig.(\ref{Fig:1 actual tumor growth}) that illustrate how the final populations of $T_S$ and $T_R$ cells change in response to variations in specific model parameters around their established baseline values. These plots offer a continuous perspective on how each parameter influences the cell populations, allowing us to observe trends and reactions over a range of adjustments. The local sensitivity analysis specifically examines the parameters $K_M$ (maximum carrying capacity) and $\kappa_1$ (rate constant) and their effects on the final populations of $T_S $ and $T_R$. According to the plots as $K_M$ increases, the final $T_{S}$ population generally increases, especially for changes above the baseline ($0\%$ change) but below the baseline, the $T_{S}$ population decreases. This confirms the positive sensitivity of $T_S$ to $K_M$. How ever on top-right plot in Fig.(\ref{Fig:1 actual tumor growth}) as $K_M$ increases, the final $T_{R}$ population generally decreases, especially for changes above the baseline but below the baseline, it increases confirming the negative sensitivity of $T_R$ to $K_M$. The final TS population shows a strong positive correlation with $\kappa_1$. As $\kappa_1$ increases, the TS population increases significantly. However, the final $T_R$ population shows a slight negative correlation with $\kappa_1$, meaning an increase in $\kappa_1$ leads to a minor decrease in TR. This effect is much less pronounced compared to its effect on TS. These plots provide a more detailed understanding of the non-linear or linear relationships between the parameters and the output variables over a range of perturbations, complementing the single-point sensitivity indices previously calculated \cite{siannis2005sensitivity, Gatenby2009}. \subsection{Global Sensitivity Analysis} In this subsection we address the inherent parametric uncertainties and structural limitations present within the baseline dynamical system, we undertook an extensive and rigorous global theoretical gradient-based sensitivity analysis. This analytical method is meticulously designed to systematically evaluate the robustness of the model while taking into account a diverse array of biological conditions that may influence its performance. To facilitate this comprehensive assessment, we made significant modifications to the foundational structural assumptions underlying the model. This involved recalibrating the equilibrium parameters to align with each specific scenario we considered, ensuring that our analysis was both relevant and contextually appropriate for the range of biological conditions being examined \cite{giudicianni2024variance}. \\ The simultaneous evaluation provides deeper insights into how variations in one parameter may affect others, ultimately contributing to a more robust understanding of the system's overall behavior and stability. In our approach, we leveraged the mathematical framework established by Mohammad et al. \cite{mohammad2021mathematical}. This framework allows for a holistic evaluation of the entire parameter manifold concurrently, presenting a distinct advantage over traditional methods. Instead of isolating parameters and analyzing them independently, our method enables us to explore their interactions and interdependencies within the broader context of the dynamical system. The bar charts in Fig.(\ref{Fig:1 global sensitive analysis}), which showcase a global sensitivity analysis partitioning the variance in the model output into contributions from individual inputs and their combinations \cite{giudicianni2024variance}. The total-order Sobol index (ST) serves as a powerful metric that quantifies the comprehensive contribution of each parameter to the variability observed in the output of a model. \begin{center} \begin{figure}[hbt!] \centering %\fbox{% %\tmpframe{ \includegraphics[width =5.7in]{image0.7.png} \includegraphics[width =5.7in]{image0.8.png} %} %} \caption{ \textit{The variance-based global sensitivity analysis (GSA) approach (Sobol) for \( T_{S} \) and \( T_{R} \) utilizes bar charts to present the first-order indices \( S1 \) and total-order indices \( ST \).}}. \label{Fig:1 global sensitive analysis} \end{figure} \end{center} It captures not only the direct effects of individual parameters but also the complex interactions that may occur with multiple other parameters \cite{Eftimie2021}. This index is critical to understanding how different factors collaboratively shape the behavior of the model’s results. In contrast, the first-order index (S1) narrows the focus to evaluate the variance that can be explained by changing only one parameter while keeping all the others constant. This means that S1 reflects a more isolated view of the influence of the parameters. In contrast, the total-order index (ST) provides a broader perspective by accounting for the variance that is not only attributable to that parameter itself, but also all its possible interactions with the entire set of parameters in the model. An elevated total-order index suggests that a given parameter plays a significant role in influencing model outputs, either on its own or through its interactions with other variables \cite{mohammad2021mathematical}. \begin{center} \begin{figure}[hbt!] \centering %\fbox{% %\tmpframe{ \includegraphics[width =3.4in]{image0.71.png} \includegraphics[width =3.0in]{image0.72.png} %} %} \caption{ \textit{The change in sensitivity for each parameter for reduced and original model with positive value indicating increased sensitivity in the reduced model, while a negative value indicates decreased sensitivity.}}. \label{Fig:1 global sensitive analysis1} \end{figure} \end{center} Parameters that significantly influence a model are essential components because they directly affect the model's behavior and the outcomes it produces. These influential parameters often have a high total-order index, which indicates that changes in their values lead to considerable adjustments in the model outputs. In contrast, parameters with a low total-order index have minimal or negligible effects on the outputs, whether their influence is direct or occurs through complex interactions with other parameters. To illustrate the sensitivity changes for each parameter clearly, we compute the differences and present them in Fig.(\ref{Fig:1 global sensitive analysis1}). A positive value here indicates increased sensitivity in the reduced model, while a negative value signifies decreased sensitivity. The Total-order Sobol index (ST) quantifies the overall effect of each parameter on output variance, encompassing both its first-order effect and all interaction effects with other parameters. By comparing the ST indices of the original and reduced models, a positive difference suggests that a parameter's overall influence on model output has increased in the reduced model. Conversely, a negative difference highlights a decrease in influence. Notably, a significant negative percentage change implies a marked reduction in the parameter's direct contribution to output variance in the reduced model compared to the original one. \\ It is important to note that, in some cases, the 'Reduced S1' values can even be negative. This may occur due to strong interactions or numerical noise, particularly when overall sensitivity is low. Nevertheless, the key takeaway is the notable reduction in the first-order impact of those parameters. However, for these specific parameters, the simplification of treating $C$ as a constant $( C^*)$ in the reduced model has drastically diminished their individual sensitivity. They are no longer as influential on the total tumor cell population on their own as they were in the full model. This could be because their direct effect was heavily mediated by $C$ in the original model, or because their influence is now primarily captured through higher-order interaction effects that became more prominent in the absence of $C$ as a dynamic variable see Appendix.(\ref{app:extra3}) for more detailed information on parameter percentage change. \section{Results and Discussions} The primary objective of this study was to emphasize the critical need for a thorough mechanistic understanding of tumor progression to inform the development of effective long-term interventions for lung cancer and to mitigate the challenge of acquired drug resistance. Through the formulation of a comprehensive quantitative framework, our research provides a robust mechanistic rationale for the frequent clinical failures associated with the use of single-agent chemotherapies and monogenic immunotherapies \cite{zhang2025mathematical}. Our findings reveal that a significant factor in these failures is the complex microenvironmental crosstalk that occurs, particularly between the immunosuppressive polarization of $M_2$ macrophages and the inactivation of effector immune cells. This interplay creates a protective sanctuary within the tumor ecosystem, allowing resistant tumor sublineages, designated as ($T_R$), to evade both the cytotoxic effects of chemotherapeutic agents and the surveillance capabilities of the host's immune system \cite{kumagai2025immunogenomic, xiong2025mathematical}.\\ To effectively reverse this established resistance, our results suggest the necessity of employing multi-target treatment strategies. Such approaches may include the integration of conventional chemotherapy ($C$) with therapies designed to induce repolarization of $M_2$ macrophages to the more tumoricidal $M_1$ phenotype. These strategies aim to enhance and preserve the functionality of effector lymphocytes ($E$) that are critical in mounting an effective anti-tumor immune response \cite{kumagai2025immunogenomic}. Additionally, by meticulously mapping the complex and non-linear feedback mechanisms inherent in tumor-immune dynamics, our model serves as a predictive tool designed to identify pivotal therapeutic tipping points. This capability allows for the optimization of combination treatment schedules and the strategic inhibition of key evolutionary pathways that underlie the dynamics of lung cancer progression \cite{xiong2025mathematical}. From a methodological perspective, this study introduces two noteworthy advancements. First, the construction of a virtual cohort surpasses conventional deterministic modeling frameworks by enabling a nuanced, population-level evaluation of the heterogeneity that exists within tumor-immune interactions \cite{xiong2025mathematical}. This multi-parameter sampling approach sheds light on how varying parameter landscapes can dictate divergent clinical trajectories, thereby enhancing our understanding of treatment responses. Second, the incorporation of optimal control theory establishes a rigorous mathematical foundation that facilitates the computation and benchmarking of the theoretical limits of treatment efficacy. This mathematical rigor ensures that our model can reliably guide future therapeutic strategies aimed at maximizing patient outcomes in lung cancer treatment \cite{kumagai2025immunogenomic}. \subsection{Optimization of Drug Resistance} In this subsection, we delve into the intricate development of an optimal control strategy aimed at addressing the mathematical modeling of cancer cell populations, with a specific emphasis on both sensitive and resistant subtypes of tumor cells. Our primary goal lies in identifying and manipulating various control variables, such as drug dosing schedules and infusion rates, that can effectively minimize the proliferation and transmission dynamics of malignant cells throughout the body. To achieve this, we formulate an optimal control problem grounded in Pontryagin's Minimum Principle (PMP) \cite{pacheco2017cancer}, which serves as a framework for determining the best possible strategies for chemotherapy administration. We introduce the control variable $ u(t) = u_0(t) $, which denotes the time-dependent infusion rate of chemotherapy drugs. This variable reflects the intensity and timing of drug delivery over the treatment course and is subject to essential physiological constraints, specifically $ 0 \leq u(t) \leq u_{\max} $. Here, $ u_{\max} $ represents the maximum feasible infusion rate that the patient's body can tolerate without significant adverse reactions. The overarching objective of our optimization is to minimize the total tumor burden, represented as the aggregate of sensitive tumor cell populations $T_S$ and resistant tumor cell populations $T_R$. In addition to this primary goal, we also place a critical focus on managing and limiting the potential toxicity associated with the administration of chemotherapeutic agents to the surrounding healthy tissues throughout the treatment duration $ [0, T_f] $. This multifaceted approach is intended to ensure a seamless balance between maximizing the therapeutic efficacy of the treatment thereby effectively reducing tumor growth and spread and minimizing the undesirable side effects and toxicity that could result from drug administration. Through the careful design of this optimal control strategy, we aim to not only enhance the effectiveness of chemotherapy in tackling tumor cells but also to contribute to the establishment of safer treatment protocols that prioritize patient well-being throughout their cancer therapy journey. In doing so, our strategy aspires to pave the way for more personalized and adaptive chemotherapy regimens that can be fine-tuned based on individual patient responses and tumor dynamics. \[J(u) = \int_{0}^{T_f} \left( A_1 T_S(t) + A_2 T_R(t) + \frac{B}{2} u^2(t) \right) dt\] where $A_1, A_2 > 0$ weight the burden of sensitive and resistant cells ($A_2 > A_1$ prioritizes suppressing resistance), and $B > 0$ penalizes high drug dosages \cite{galete2025mathematical}. Let $x = (T_S, T_R, E, M_1, M_2, C)^T$ be the state vector and $\lambda = (\lambda_1, \lambda_2, \lambda_3, \lambda_4, \lambda_5, \lambda_6)^T$ be the adjoint (costate) vector. The Hamiltonian $H$ is defined as: \begin{eqnarray} H(x, u, \lambda, t) = A_1 T_S + A_2 T_R + \frac{B}{2} u^2 + \lambda_1 \frac{dT_S}{dt} + \lambda_2 \frac{dT_R}{dt} + \lambda_3 \frac{dE}{dt} + \lambda_4 \frac{dM_1}{dt} + \lambda_5 \frac{dM_2}{dt} + \lambda_6 \frac{dC}{dt}\label{eq1.12} \end{eqnarray} By Pontryagin's Minimum Principle, the adjoint variables satisfy $\frac{d\lambda_i}{dt} = -\frac{\partial H}{\partial x_i}$ with transversality conditions $\lambda_i(T_f) = 0$ for $i = 1, \dots, 6$: \begin{eqnarray} \begin{aligned} \frac{d\lambda_1}{dt} =& -A_1 - \lambda_1 \left[ r_S \left(1 - \frac{2T_S + T_R}{K_T}\right) + \alpha M_2 - \delta_1 M_1 - \kappa_1 E - \psi C - \mu(M_2, C) \right] \\ & - \lambda_2 \left[ -\frac{r_R T_R}{K_T} + \mu(M_2, C) \right] - \lambda_3 \frac{\rho_0 a_0 E}{(a_0 + T_S + T_R)^2} + \lambda_4 \delta M_1 - \lambda_5 \frac{\rho_2 a_2}{(a_2 + T_S + T_R)^2}, \\ \frac{d\lambda_2}{dt} =& -A_2 - \lambda_1 \left( \eta - \frac{r_S T_S}{K_T} \right) - \lambda_2 \left[ r_R \left(1 - \frac{T_S + 2T_R}{K_T}\right) + \alpha M_2 - \delta_2 M_1 - \kappa_2 E - \eta - \varepsilon \psi C \right] \\ & - \lambda_3 \frac{\rho_0 a_0 E}{(a_0 + T_S + T_R)^2} + \lambda_4 \delta M_1 - \lambda_5 \frac{\rho_2 a_2}{(a_2 + T_S + T_R)^2}, \\ \frac{d\lambda_3}{dt} =& \lambda_1 \kappa_1 T_S + \lambda_2 \kappa_2 T_R - \lambda_3 \left[ \frac{\rho_0(T_S + T_R)}{a_0 + T_S + T_R} - \mu_E - \phi M_2 - 2\theta E \right] - (\lambda_4 - \lambda_5) \frac{\rho_1 a_1 M_2}{(a_1 + E)^2}, \\ \frac{d\lambda_4}{dt} =& \lambda_1 \delta_1 T_S + \lambda_2 \delta_2 T_R - \lambda_4 \left[ p_{M1} \left(1 - \frac{2M_1 + M_2}{K_M}\right) - \mu_{M_1} - \delta(T_S + T_R) \right] + \lambda_5 \frac{p_{M2} M_2}{K_M}, \\ \frac{d\lambda_5}{dt} =& -\lambda_1 (\alpha - \mu_1) T_S - \lambda_2 (\alpha T_R + \mu_1 T_S) + \lambda_3 \phi E + \lambda_4 \frac{p_{M1} M_1}{K_M} - (\lambda_4 - \lambda_5) \frac{\rho_1 E}{a_1 + E} \\ & - \lambda_5 \left[ p_{M2} \left(1 - \frac{M_1 + 2M_2}{K_M}\right) - \mu_{M_2} \right], \\ \frac{d\lambda_6}{dt} =& \lambda_1 (\psi + \mu_2) T_S + \lambda_2 (\varepsilon \psi T_R - \mu_2 T_S) + \lambda_6 k_C. \label{eq1.1} \end{aligned} \end{eqnarray} To determine the optimal dosing schedule $u^*(t)$, we apply the optimality condition $\frac{\partial H}{\partial u} = 0$: \[\frac{\partial H}{\partial u} = B u(t) + \lambda_6(t) = 0 \implies u(t) = -\frac{\lambda_6(t)}{B}\] Taking into account the upper and lower bounds $0 \le u(t) \le u_{\max}$, the optimal drug administration rate is given by: \[u^*(t) = \min \left( u_{\max}, \max \left( 0, -\frac{\lambda_6(t)}{B} \right) \right)\] This forward-backward sweep system, which involves integrating state differential equations forward from the initial time $t=0$ and costate equations backward from the final time $t=T_f$, offers a precise adaptive protocol designed to minimize the development of drug resistance. This resistance is primarily driven by the activity of $M_2$ macrophages and the dose-dependent mutation rate represented as $\mu(M_2, C)$. \begin{center} \begin{figure}[hbt!] \centering %\fbox{% %\tmpframe{ \includegraphics[width = 6.3in]{image0.9.png} %\includegraphics[width = 6.34in]{image0.8.png} %} %} \caption{ \textit{the Forward-Backward Sweep Method (FBSM) accelerated with 4th-Order Runge-Kutta (RK4) integration and a convex convergence update to solve the optimal control system.}}. \label{Fig:1 optimal control system} \end{figure} \end{center} As highlighted by Galete et al. \cite{galete2025mathematical}, a thorough assessment of the sensitivity and robustness of the optimized model presented in Eq.\eqref{eq1.1} requires us to delve into the parameter influence profiles illustrated in Figures \ref{Fig:1 optimal control system} and \ref{Fig:1 optimal control system 1}. These sensitivity plots vividly demonstrate how critical biological and kinetic parameters impact the optimal control trajectory $u^*(t)$ and the effectiveness of strategies aimed at mitigating the emergence of drug-resistant tumor sub-populations. Notably, Fig.(\ref{Fig:1 optimal control system 1}) provides detailed insights into the optimal control profiles corresponding to the parameter $w_3$. In panel A, the graph showcases the optimal chemotherapy infusion rate, denoted as $u^{*}(t)$, plotted over time for a range of different $w_3$ values. \begin{center} \begin{figure}[hbt!] \centering %\fbox{% %\tmpframe{ \includegraphics[width = 6.3in]{image1.1.png} %\includegraphics[width = 6.34in]{image0.8.png} %} %} \caption{ \textit{Varying the control cost weight $w_3$ directly modulates the trade-off between drug toxicity tolerance and tumor suppression aggressiveness: lower values of $w_3$ permit aggressive, near-maximal dosing, while higher values force conservative schedules to minimize systemic toxicity.}}. \label{Fig:1 optimal control system 1} \end{figure} \end{center} This visualization allows for a clearer understanding of how variations in \(w_3\) influence the chemotherapy protocols over the treatment period. The $w_3$ parameter dictates the penalty associated with using chemotherapy. As $w_3$ increases (meaning a higher cost or penalty for the drug), the optimal infusion rate generally decreases, indicating a more conservative dosing strategy to minimize the control cost. This plot shows the combined population of sensitive ($T_{S}$) and resistant ($T_{R}$) tumor cells over time for each $w_{3}$ value. It directly reflects the effectiveness of the different chemotherapy dosing regimens. Lower $w_{3}$ values, which typically lead to higher drug doses, often result in a more significant reduction in the total tumor burden. Panel C illustrates the concentration of chemotherapy ($C(t)$) in the system over time, corresponding to the different $w_{3}$ values. The drug concentration is a direct consequence of the applied dosing strategy (from Panel A) and the drug's elimination rate. Higher doses generally lead to higher peak concentrations and longer exposure times. However panel D shows the relationship between the control weight $w_3$ (on a logarithmic scale) and the total cumulative dose of chemotherapy administered over the entire treatment period. It clearly demonstrates the trade-off: as the penalty for using chemotherapy ($w3$) increases, the total amount of drug administered decreases, aiming to find a balance between tumor control and treatment cost or toxicity. \begin{theorem}{(Characterization of Optimal Control)} An optimal control $u(t)$ and corresponding state solutions $(T_S^*, T_R^*, E^*, M_1^*, M_2^*, C^*)$ that minimize $J(u)$ over $U_{ad}$ exist. Moreover, there exist adjoint functions $\lambda_i(t)$ ($i=1,\dots,6$) satisfying the adjoint system given in Eq.\eqref{eq1.1} along with the transversality conditions: \[\lambda_i(T_f) = 0 \quad \text{for } i = 1, 2, 3, 4, 5, 6\] Furthermore, the optimal control $u(t)$ is characterized by: \[u^*(t) = \min \left( u_{max}, \max \left( u_{min}, -\frac{\lambda_6(t)}{B} \right) \right)\] \end{theorem} \begin{proof} To prove the Optimal Control Theorem, we rigorously apply Pontryagin's Minimum Principle (PMP), which serves as a foundational tool in optimal control theory. We begin by establishing the existence of an optimal control $u^*$ within the set of admissible controls $U_{ad}$, as well as the associated state trajectories that arise from this control. This verification is grounded in the well-documented existence results provided by Fleming and Rishel in their seminal work from 1975. It is essential to note that the following conditions must be met in order to ensure the applicability of these principles: \begin{itemize} \item Non-emptiness of $U_{ad}$ with the control set $U_{ad} = \{ u(t) \in L^\infty(0, T_f] \mid u_{min} \le u(t) \le u_{max} \}$ is closed, bounded, and convex in $L^\infty(0, T_f]$. \item Stating that the right-hand side of state system \eqref{eq:model} is continuously differentiable ($\mathcal{C}^1$) with respect to state variables and linear in $u(t)$. Thus, state solutions remain bounded on $[0, T_f]$. \item Conversely, the integrand of the objective functional $f(t, \mathbf{x}, u) = A_1 T_S + A_2 T_R + \frac{1}{2} B u^2$ is strictly convex with respect to $u(t)$ on $U_{ad}$. \item Then for the convex lower bound, the integrand satisfies $f(t, \mathbf{x}, u) \ge \frac{1}{2} B u^2$ for $B > 0$, ensuring the quadratic penalty dominates as $u$ varies. \end{itemize} Hence, an optimal control pair $(u^*, \mathbf{x}^*)$ exists. Formulating the Hamiltonian, we now let the variables $\mathbf{x} = (T_S, T_R, E, M_1, M_2, C)^T$ denote the vector of state variables, and $\boldsymbol{\lambda} = (\lambda_1, \lambda_2, \lambda_3, \lambda_4, \lambda_5, \lambda_6)^T$ denote the vector of adjoint (costate) variables, then the Hamiltonian $H(t, \mathbf{x}, u, \boldsymbol{\lambda})$ is defined by Eq.(\ref{eq1.12}). To derive the adjoint equations in Eq.(\ref{eq1.1}) along with the transversality conditions, we utilize Pontryagin's Minimum Principle. According to this principle, the adjoint equations are determined by the negative partial derivatives of the Hamiltonian $ H $ with respect to each corresponding state variable, evaluated at the optimal states. \[\frac{d\lambda_i}{dt} = -\frac{\partial H}{\partial x_i}, \quad i = 1, \dots, 6\] In order to develop a control characterization aimed at minimizing the Hamiltonian $H$ with respect to the control variable $u(t)$ within the admissible set $U_{ad}$, we first calculate the partial derivative of the Hamiltonian. This step is crucial for understanding how changes in the control input affect the overall system dynamics and performance measures. By analyzing this derivative, we can identify optimal control strategies that lead to a reduced value of $ H$. \[\frac{\partial H}{\partial u} = B u(t) + \lambda_6(t)\] Since \( B > 0 \), the second derivative \( \frac{\partial^2 H}{\partial u^2} = B > 0 \) confirms that \( H \) is strictly convex with respect to \( u \). Evaluating at the interior point where \( \frac{\partial H}{\partial u} = 0 \): \[B u(t) + \lambda_6(t) = 0 \implies u(t) = -\frac{\lambda_6(t)}{B}\] \[u^*(t) = \min \left( u_{max}, \max \left( u_{min}, -\frac{\lambda_6(t)}{B} \right) \right)\] This completes the proof. \end{proof} \subsubsection{Existence of Optimal Control} In this subsection, we aim to establish the existence of an optimal control solution for the problem defined in Eq.(\ref{eq1.1}). To effectively achieve this goal, we will first introduce a fundamental theorem that lays out the essential theoretical framework and the specific conditions that need to be satisfied for such a solution to exist. This theorem will be pivotal in our analysis, acting as a crucial stepping stone throughout the proof process. It will help us systematically navigate the complexities of the optimization problem, ensuring that we consider all relevant factors and nuances that may influence the outcome of our control solution. \begin{theorem} Let the admissible control set $\mathcal{U}$ be compact and convex in $L^\infty([0, T_f], \mathbb{R}^m)$. Given non-negative initial conditions for state variables $(T_S, T_R, E, M_1, M_2, C)$, there exists an optimal control vector $u(t) \in \mathcal{U}$ and an associated state vector trajectory $(T_S^*, T_R^*, E^*, M_1^*, M_2^*, C^*)$ on $t \in [0, T_f]$ that satisfies: \[J(u^*(t)) = \inf_{u \in \mathcal{U}} J(u(t))\] \end{theorem} \begin{proof} To prove the existence of an optimal control $u^*(t) \in \mathcal{U}$, we apply the Filippov-Cesari Existence Theorem. Firstly, we verify the compactness and convexity of the control space where the admissible control set is defined by \[\mathcal{U} = \left\{ u \in L^\infty([0, T_f], \mathbb{R}^m) \; \Big\vert{} \; 0 \le u_j(t) \le u_j^{\max}, \, \forall t \in [0, T_f] \right\}\] which is closed, bounded, and convex in $L^\infty([0, T_f], \mathbb{R}^m)$, and thus weakly compact in $L^2([0, T_f], \mathbb{R}^m)$.\\ Then for the positivity and invariance of states, the vector field $f(t, x, u)$ satisfies $\left. f_i(t, x, u) \right\vert{}_{x_i = 0} \ge 0$ for all state components $x = (T_S, T_R, E, M_1, M_2, C)$. According to Nagumo's Theorem, the non-negative orthant $\mathbb{R}^6_+$ is positively invariant.\\ Since cell dynamics is bounded above by logistic carrying capacities $K$ and natural or clearance death rates, all state variables remain uniformly bounded in a compact region $\Omega \subset \mathbb{R}^6_+$ for all $t \in [0, T_f]$.\\ The right-hand side of the state ODE system is linear with respect to the control vector $u(t)$, taking the form: \[f(t, x, u) = g(t, x) + \sum_{i=1}^{m} h_i(t, x) u_i(t)\] However, using the objective functional integrand as follows from Lassounon et al. \cite{lassounon2026optimal}, then $L(t, x, u) = \sum_{k} w_k x_k + \sum_{j} \frac{B_j}{2} u_j^2$ is continuous and strictly convex with respect to $u$ due to the quadratic control terms ($\frac{B_j}{2} u_j^2$). In full for \( j=1,2,3 \), it is imperative that the approach adheres to several key conditions. Firstly, prevention controls must be strategically designed to significantly reduce the transmission rate of the disease. Secondly, cancer treatment should prioritize timely and effective interventions, specifically focusing on the administration of second-line drugs that have proven effective in managing the disease. Furthermore, the control of cancer treatment must specifically target the most resistant forms of the disease, which present formidable challenges due to their high mortality rates and the limited availability of viable treatment options.\\ According to various studies \cite{lassounon2026optimal, okosun2013impact}, it becomes evident that a comprehensive and integrated strategy is essential for both prevention and treatment, ultimately enhancing the overall efficacy of the healthcare response. Additionally, continuous research and the refinement of treatment protocols are crucial to adapt to the changing landscape of cancer challenges, ensuring that healthcare responses remain effective and relevant. Then from the conditions we mentioned the theorem is expanded to give: \[L(t,x,u)\geq \min \left(\frac{w_{1}}{2},\frac{w_{2}}{2}, \frac{w_{3}}{2}\right)(u_{1}^{2}+u_{2}^{2}+u_{3}^{2})-\frac{w_{2}}{2}\] \[\Longrightarrow L(t,x,u)\geq \min \left(\frac{w_{1}}{2},\frac{w_{2}}{2}, \frac{w_{3}}{2}\right)\Vert u_{1}^{2}+u_{2}^{2}+u_{3}^{2}\Vert ^{2}-\frac{w_{2}}{2}\] The function $L(t,x,u)$ is bounded if and only if the parameter $\zeta$ satisfies $\zeta = \frac{1}{2} \min \{w_1, w_2, w_3\}$. Since the objective functional $J(u)$ is bounded from below on the admissible control space $\mathcal{U}$, there exists a minimizing sequence $\{u^{(k)}\} \subset \mathcal{U}$ satisfying \[\lim_{k \to \infty} J(u^{(k)}) = \inf_{u \in \mathcal{U}} J(u)\] By weak compactness, $\{u^{(k)}\}$ has a subsequence (still $\{u^{(k)}\}$) that converges weakly in $L^2([0, T_f])$ to some $u^* \in \mathcal{U}$.The corresponding state solutions $\{x^{(k)}\}$ are uniformly bounded and equicontinuous. By the Arzelà-Ascoli Theorem, $x^{(k)} \to x^*$ uniformly on $[0, T_f]$.Since $L(t, x, u)$ is convex with respect to $u$, $J(u)$ is weakly lower semicontinuous. Therefore: \[J(u^*) \le \liminf_{k \to \infty} J(u^{(k)}) = \inf_{u \in \mathcal{U}} J(u)\] Thus, $J(u^*) = \min_{u \in \mathcal{U}} J(u)$, confirming that $u^*(t)$ is an optimal control with corresponding optimal state trajectory $x^*(t) = (T_S^*, T_R^*, E^*, M_1^*, M_2^*, C^*)$ \end{proof} \subsection{Treatment and Resistance Scenario} The treatment of lung cancer encompasses a variety of approaches, including chemotherapy, targeted therapy, and immunotherapy. Each of these treatment modalities aims to combat the disease through different mechanisms, yet a significant obstacle that arises in managing lung cancer is the development of drug resistance. This phenomenon can be categorized into two primary types, which includes intrinsic resistance, which is present from the outset of treatment, and acquired resistance, which develops over time as the cancer evolves \cite{emran2022multidrug}. \\ The mechanisms behind this resistance are multifaceted and complex. For instance, drug efflux, a process where cancer cells actively expel therapeutic agents, can diminish the effectiveness of treatments. Additionally, mutations may occur in the genes targeted by these drugs, altering their effectiveness. Moreover, changes within the tumor microenvironment itself, such as an altered supply of nutrients or oxygen, can create a setting that further promotes resistance. Understanding these mechanisms is crucial for developing strategies to overcome resistance and improve treatment outcomes in lung cancer patients \cite{lei2023understanding, khalaf2021aspects}. To overcome resistance in cancer treatment, several strategies have been proposed. These include using drugs that operate through different molecular mechanisms to enhance tumor cell killing while minimizing the likelihood of resistance. Additionally, utilizing checkpoint inhibitors or monoclonal antibodies can bolster the immune response against cancer cells. These approaches aim to optimize the effectiveness of cancer treatments and improve clinical outcomes by addressing the complex mechanisms of drug resistance \cite{li2025drug}. To evaluate therapeutic impact, the systemic or intratumoral drug concentration $C(t)$ defined in Eq.(\ref{eq:C}) must first be calculated using an integrating factor $e^{k_C t}$, the differential equation is solved given the initial condition $C(0) = C_0$ to get: \[C(t) = C_0 e^{-k_C t} + \frac{u_0}{k_C} \left(1 - e^{-k_C t}\right)\] Steady state concentration ($C_{ss}$) is an important pharmacokinetic concept that is reached as time approaches infinity ($t \to \infty$). At this point, the rate of drug elimination balances the rate of drug input, leading to a stable concentration of the drug in the system. Understanding $C_{ss}$ is essential for developing effective dosing strategies, ensuring that therapeutic levels are achieved while minimizing the risk of toxicity. According to Pradhan et al. \cite{pradhan1995steady} the steady state of drug of concentration is defined from: \[C_{ss} = \lim_{t \to \infty} C(t) = \frac{u_0}{k_C}\] Introducing drug Clearance and elimination half-life ($t_{1/2}$), which defines the time required for drug concentration to decay by $50\%$ in the absence of input ($u_0 = 0$). The concept of clearance refers to the volume of plasma from which the drug is completely removed per unit time. Understanding these parameters is crucial for determining appropriate dosing regimens and assessing the pharmacokinetics of a drug within the body. To calculate the elimination half-life, the formula is: \[ t_{1/2} = \frac{0.693 \cdot V_d}{Cl}=\frac{ln(2)}{k_{C}} \] In pharmacokinetics, \(V_d\) refers to the volume of distribution, which indicates how extensively a drug disperses throughout the body relative to the concentration of the drug in the blood plasma. Clearance (\(Cl\)), on the other hand, is a measure of the body's efficiency in eliminating the drug, accounting for the rate at which the drug is removed from systemic circulation. Adjusting either the clearance or the volume of distribution can have a profound effect on the half-life of a medication, a critical factor in the therapeutic management of patients. The half-life is the time it takes for the plasma concentration of a drug to reduce by half, and this metric is crucial for determining dosing intervals and durations of treatment to maximize efficacy while minimizing potential side effects \cite{pradhan1995steady, li2025drug}. Utilizing the pharmacokinetic equations that detail the relationships among these variables, we can perform an in-depth analysis of drug concentration levels in patients. This analysis is particularly important for detecting strains of pathogens that have developed resistance to standard treatment protocols. Such comprehensive evaluations are essential not only for identifying the existence of drug-resistant strains but also for understanding how varying infusion rates influence the overall concentration of the drug circulating in the bloodstream \cite{fernandez2019optimal}. This knowledge allows clinicians to tailor infusion rates and dosing schedules to achieve optimal therapeutic levels of the drug. Moreover, grasping these pharmacokinetic dynamics is vital for effective therapeutic drug monitoring. By continuously assessing drug levels in the body, healthcare providers can adjust treatment plans in real-time, ensuring that drug concentrations remain within the therapeutic window. This personalized approach improves patient care and enhances treatment efficacy, ultimately contributing to better health outcomes \cite{ramos2021battling}. \begin{center} \begin{figure}[hbt!] \centering %\fbox{% %\tmpframe{ \includegraphics[width =3.3in]{image4.2.png} \includegraphics[width =3.3in]{image4.21.png} \includegraphics[width =3.3in]{image4.22.png} \includegraphics[width =3.3in]{image4.23.png} %} %} \caption{\textit{The relationship between steady-state concentration and drug concentration over time, along with the impact of the infusion rate \( u_0 \) on the drug concentration curve and steady-state levels, as well as the time needed to achieve these concentrations.}}. \label{Fig:1 Drug concentration} \end{figure} \end{center} To thoroughly analyze systemic drug accumulation dynamics, we created detailed numerical profiles that illustrate the pharmacokinetic response across a range of infusion rates as shown in Fig.(\ref{Fig:1 Drug concentration}). Consistent with our analytical model, we observed that the steady-state concentration ($C_{ss}$) maintains a strictly proportional linear relationship with the infusion rate of zero-order ($u_0$). For example, when the infusion rate doubles from $100$ to $200$, there is a corresponding increase in steady-state concentration, which also doubles to $C_{ss} = 2000.00$. Importantly, adjusting the infusion rate $u_0$ influences only the amplitude of the equilibrium concentration; it does not change the characteristic time required to reach steady-state. Our analysis revealed that the system achieves $90\%$ of the steady-state concentration at approximately $t \approx 23.03$ hours. This finding underscores that the rate at which the drug accumulates to reach steady-state is determined solely by the elimination rate constant ($k_C$) and is not affected by the infusion rate. Understanding and quantifying these steady-state percentages, along with the associated clearance timescales, is crucial to the development of optimal dosing regimens. This knowledge helps prevent sublethal selection windows in treatment, ensuring that effective washout periods are attained and maintaining the overall efficacy of the therapeutic regimen \cite{pradhan1995steady}. The steady-state concentration (Css) is reached when the rate of drug administration equals the rate of elimination, producing a relatively constant plasma drug level. Maintaining high fractions of Css ($>80\%$) ensures maximum target engagement, particularly relevant for tyrosine kinase inhibitors, immune checkpoint inhibitors, and conventional chemotherapeutic used in lung cancer. High steady-state drug concentrations ($>80\%$ Css) can maximize initial target engagement but also accelerate the selection and emergence of drug-resistant, persistent, or dormant tumor cells \cite{garg2024emerging, gottesman2002mechanisms}. \subsection{Model Comparison For Standard Treatment} According to Injerd et al. \cite{injerd2017mathematical}, Akaike’s Information Criterion (AIC) was utilized to systematically compare and evaluate the relative quality of different models for each individual patient. This criterion serves to optimize the trade-off between statistical goodness of fit and model complexity, ensuring the selection of a parsimonious model that adequately explains the patient data while minimizing over-parameterization. The AIC can be written as \begin{eqnarray} AIC=T\left(\ln{(2\pi)}+\ln{(SSE/T)}+1\right)+2p \end{eqnarray} where $T$ represents the total number of observations, $SSE$ is the sum of squared errors measuring how far the model's predictions are from the actual values and $p$ is the total number of estimated parameters in the model. When dealing with small samples, as is often the case in oncological modeling, AIC may lead to over fitting. Although standard AIC performs reliably with large sample sizes, it exhibits a known negative bias in small sample regimes such as those frequently encountered in oncological and clinical modeling, leading to a tendency to select over fitted and overly complex models. To mitigate this issue, we adopt the small sample corrected Akaike Information Criterion ($\text{AIC}_c$) proposed by Hurvich and Tsai (1989) \cite{hurvich1989regression}. To prevent overfitting in finite-sample regimes, $\text{AIC}_c$ modifies the standard $\text{AIC}$ metric with a correction factor inversely proportional to the residual degrees of freedom, $T - p - 1$: \[\text{AIC}_c = \text{AIC} + \frac{2p(p + 1)}{T - p - 1}\] Substituting the explicit formula for AIC yields the following: \[\text{AIC}_c = T \ln\left(\frac{\text{SSE}}{T}\right) + 2p + \frac{2p(p + 1)}{T - p - 1}\] As $T \to \infty$, the correction term vanishes, causing $\text{AIC}_c$ to converge asymptotically to a standard AIC. Specifically, when the sample size $T$ is small relative to the number of parameters $p$ conventionally defined as a ratio where $T / p < 40$ for standard AIC under-penalizes parameter proliferation \cite{giudicianni2024variance}. In these data-limited regimes, $\text{AIC}_c$ applies a progressively steeper penalty as $p$ approaches $T-1$, preventing the inclusion of spurious parameters and safeguarding against overfitting in clinical modeling \cite{hurvich1989regression}.To complement $\text{AIC}_c$ and eliminate selection bias during model comparison, we additionally employ the Bayesian Information Criterion (BIC), formulated as: \[\text{BIC} = T \ln\left(\frac{\text{SSE}}{T}\right) + p \ln(T)\] While $\text{AIC}_c$ seeks to minimize relative information loss (Kullback-Leibler divergence) to maximize predictive accuracy, BIC estimates the posterior probability of a model being the true data-generating process \cite{neath2012bayesian}. Because BIC scales its penalty term by $\ln(T)$ rather than a constant factor of $2$, it imposes a heavier penalty on model complexity for any sample size $T \ge 8$ (since $\ln(8) \approx 2.08 > 2$). The criteria for both model from Eq.(\ref{eq:TS}) to Eq.(\ref{eq:C}) and the reduced model from Eq.(\ref{eq:RTS}) to Eq.(\ref{eq:RM2}) against the same synthetic noisy data derived from first model's simulation is calculated and gives the results shown in below Table.(\ref{tab:AIC}). \begin{table}[htbp] \centering \caption{Model Comparison} \label{tab:AIC} \footnotesize \renewcommand{\arraystretch}{1.2} \begin{tabular}{lllll} \hline Criterion & Model $1$ Value & Reduced Model & Preferred \\ & & Value & Model \\ \hline AIC & $81778.02$ & $78900.45$ & Reduced\\ AICc& $81778.73$ & $78901.26$ & Reduced\\ BIC & $81970.23$ & $79081.00$ & Reduced\\ \hline \end{tabular} \end{table} The AIC criterion suggests that the Reduced Model is a better fit for the (noisy) data, implying that the simplification of the drug concentration dynamics ($C$) provides a more efficient description of the system without losing significant explanatory power. We have performed a comprehensive comparison between Model $1$ and the Reduced Model using three key information criteria AIC, AICc, and BIC. The results consistently point to the same conclusion, all three information criteria consistently favor the reduced Model over model $1$. In practical terms, this implies that for the given set of parameters and conditions, the dynamic evolution of the drug concentration $C$ (as modeled in Model 1) d, given the set of parameters and conditions, the dynamic evolution of the drug concentration (as modeled in Model 1) does not provide enough additional explanatory power to justify the increased complexity of the model \cite{injerd2017mathematical}. The simpler reduced Model provides a more efficient and statistically sound description of the system's dynamics. \section{Conclusion} In this comprehensive study, we developed and rigorously analyzed a sophisticated compartmental mathematical framework specifically designed to investigate the intricate population dynamics, evolutionary selection mechanisms, and optimal therapeutic control strategies related to drug resistance in lung cancer \cite{galete2025mathematical}. Our approach allowed us to delve deeply into the mechanisms by which lung cancer cells adapt and develop resistance against treatment, which is crucial for informing more effective therapeutic interventions. To predict tumor growth or regression analytically, we utilized the next-generation matrix approach, a powerful tool in mathematical biology that facilitates the analysis of complex systems. Through this methodology, we was able to derive the basic reproduction number ($R_0$) analytically. This parameter is essential as it represents the threshold governing cellular fitness, evolutionary persistence, and the overall growth dynamics of tumors \cite{galete2025mathematical, bracci2014immune}.\\ The significance of $R_0$ lies in its ability to act as a predictor for whether a population of cancer cells will decline or proliferate, thereby guiding therapeutic decisions. In our analysis, we also conducted both local and global asymptotic stability analyses, which were formally established in Theorem.(\ref{theorem1.3}). The findings are particularly important, as they indicate a critical threshold which suggest the treatment strategies in place are sufficiently effective to control tumor growth, and the cancer population is likely to diminish over time \cite{bracci2014immune}. Moreover, our model reveals a forward bifurcation at the threshold value of $R_0>1$. This bifurcation underscores the critical nature of sustaining $R_0$ strictly below unity, reinforcing that this condition is an indispensable prerequisite for both eradicating the disease and preventing its resurgence \cite{ahmad2024mathematical}. This insight paves the way for the development of therapeutic strategies aimed at not only controlling existing tumors but also preventing the emergence of drug-resistant cancer cell populations. Consequently, the design and implementation of clinical interventions require a strategic optimization approach to effectively target and suppress both drug-sensitive and drug-resistant sub-lineages of malignant cells. This optimization holds the potential to significantly decrease the overall burden of lung cancer on patients and healthcare systems alike. By integrating principles from optimal control theory with the dynamics of lung cancer populations, this research tackles fundamental open questions in the fields of mathematical oncology and disease modeling, offering novel insights into treatment efficacy and patient outcomes \cite{touray2025combating}. The framework developed in this study meticulously incorporates a range of crucial real-world clinical variables. These include treatment adherence rates, which reflect how well patients are following prescribed therapeutic regimens; age-structured host vulnerabilities, acknowledging the varying responses of different age groups to cancer treatments; early detection protocols, which play a critical role in identifying lung cancer cases at manageable stages; and comprehensive systemic therapeutic coverage, which ensures that all potential therapeutic options are utilized effectively. By accounting for these factors, the framework addresses significant structural limitations that have historically hindered conventional cancer models \cite{ashrafi2022current}. \\ In addition to its immediate contributions to mathematical findings, this flexible framework establishes a solid theoretical foundation that promotes interdisciplinary collaboration among researchers and clinicians. It is designed to foster future translational research efforts, facilitating the exchange of knowledge and resources across national and international institutions. This collaborative approach could ultimately enhance the effectiveness of lung cancer management strategies and improve patient care on a global scale \cite{ahmad2025analytical}. Despite the valuable insights offered by this study, several methodological limitations must be taken into account. One significant challenge arises from the limited availability of longitudinal time series data, which makes it impossible to conduct direct parameterization for the distinct sub-clusters identified within the tumor population \cite{chen2023macrophages}.\\ To overcome this obstacle, parameter estimation was performed using cross-sectional tumor data while adhering to steady state assumptions. Although necessary, this approach introduces certain constraints; however, it is important to note that the model effectively captures the essential mechanisms that drive the progression of lung cancer. Moreover, it serves as a flexible foundational framework that can be refined and updated as longitudinal gene expression profiles from patients become accessible, thus enhancing the model's accuracy over time. Looking to the future, there is a critical opportunity to build on this framework, especially when considering the significant impact that targeted therapies have on modulating host immunity. By extending this model to incorporate therapeutic interactions directly within the tumor microenvironment, researchers can work towards optimizing adaptive dosing schedules \cite{young2024data, nam2024harnessing}. This enhancement could lead to more personalized treatment strategies, ultimately improving outcomes for patients with lung cancer as new insights into the dynamics between the tumor, the immune system, and therapeutic agents are integrated into the model. % Start appendix \appendix \section*{APPENDIX} \section{Abbreviations}\label{app:extra1} This is the first appendix section. %Appendices are labeled alphabetically by default. \begin{table}[htbp] \centering \caption{Abbreviations} \label{tab:app1} \footnotesize \renewcommand{\arraystretch}{0.8} \begin{tabular}{lrrr} \hline Abbreviation & Meaning~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\\ \hline AICc& Small Sample Akaike Information Criterion\\ AIC & Akaike Information Criteria~~~~~~~~~~~~~~~~~~~~~~\\ AUC & Area Under the Curve~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\\ BIC & Bayesian Information Criteria~~~~~~~~~~~~~~~~~~~\\ GSA & Global Sensitivity Analysis~~~~~~~~~~~~~~~~~~~~~~~\\ LSA & Local Sensitivity Analysis~~~~~~~~~~~~~~~~~~~~~~~~\\ S1 & First-order Sobol Index~~~~~~~~~~~~~~~~~~~~~~~~~~~\\ ST & Total-order Sobol Index~~~~~~~~~~~~~~~~~~~~~~~~~~~\\ \hline \end{tabular} \end{table} \section{Additional Data} This is the second appendix section. Appendices are labeled alphabetically by default. \subsection{Next Generation Matrix}\label{app:extra2} Calculate the reproductive number ($R_0$) of the system using the Next Generation Matrix method following the method of van den et al. \cite{van2002reproduction}, we focus on the compartments of tumor cells (infected or disease-carrying).Then $T_S$ (sensitive tumor cells) and $T_R$ (resistant tumor cells) at the tumor-free equilibrium (TFE), $T_S^* = 0$ and $T_R^* = 0$, while the non-tumor variables settle at their respective steady-state values $(E^*, M_1^*, M_2^*, C^*)$, where $\mu^* = \mu(M_2^*, C^*) = \mu_0 + \mu_1 M_2^* + \mu_2 C^*$. We start by formulating the rate vectors $\mathcal{F}$ and $\mathcal{V}$ by splitting the rate of change of the tumor compartments: \[\mathbf{x} = \begin{pmatrix} T_S \\ T_R \end{pmatrix}\] into $\mathcal{F}_i$ the rate of proliferation of new tumor cells $\mathcal{V}_i = \mathcal{V}_i^- - \mathcal{V}_i^+$: the net rate of transfer outside the compartment $i$ (clearance, death, mutation/transitions). For the system: \[\mathcal{F}(T_S, T_R) = \begin{pmatrix} r_S T_S \left(1 - \frac{T_S + T_R}{K_T}\right) + \alpha M_2 T_S \\ r_R T_R \left(1 - \frac{T_S + T_R}{K_T}\right) + \alpha M_2 T_R \end{pmatrix}\] \[\mathcal{V}(T_S, T_R) = \begin{pmatrix} (\delta_1 M_1 + \kappa_1 E + \psi C + \mu(M_2, C)) T_S - \eta T_R \\ (\delta_2 M_1 + \kappa_2 E + \varepsilon \psi C + \eta) T_R - \mu(M_2, C) T_S \end{pmatrix}\] Linearize at the Tumor-Free Equilibrium ($\hat{E}^{*}$) by evaluating the Jacobians at $T_S^* = 0, T_R^* = 0$ and the non-zero $\hat{E}^{*}$ states $(E^*, M_1^*, M_2^*, C^*)$: \begin{itemize} \item[1.] Proliferation Matrix $F$: \[F = \begin{pmatrix} \left.\frac{\partial \mathcal{F}_1}{\partial T_S}\right\vert{}_{\hat{E}^{*}} & \left.\frac{\partial \mathcal{F}_1}{\partial T_R}\right\vert{}_{\hat{E}^{*}} \\ \left.\frac{\partial \mathcal{F}_2}{\partial T_S}\right\vert{}_{\hat{E}^{*}} & \left.\frac{\partial \mathcal{F}_2}{\partial T_R}\right\vert{}_{\hat{E}^{*}} \end{pmatrix} = \begin{pmatrix} r_S + \alpha M_2^* & 0 \\ 0 & r_R + \alpha M_2^* \end{pmatrix}\] For simplicity, let $\beta_S = r_S + \alpha M_2^*$ and $\beta_R = r_R + \alpha M_2^*$. Then: \[F = \begin{pmatrix} \beta_S & 0 \\ 0 & \beta_R \end{pmatrix}\] \item[2.] Transition or clearance matrix $V$ using the defined expressions $d_S = \delta_{1} M_{1}^* + \kappa_1 E^* + \psi C^* + \mu^*$ and $d_{R} = \delta_{2} M_{1}^* + \kappa_{2} E^* + \varepsilon \psi C^* + \eta$: \[V = \begin{pmatrix} \left.\frac{\partial \mathcal{V}_1}{\partial T_S}\right\vert{}_{\hat{E}^{*}} & \left.\frac{\partial \mathcal{V}_1}{\partial T_R}\right\vert{}_{\hat{E}^{*}} \\ \left.\frac{\partial \mathcal{V}_2}{\partial T_S}\right\vert{}_{\hat{E}^{*}} & \left.\frac{\partial \mathcal{V}_2}{\partial T_R}\right\vert{}_{\hat{E}^{*}} \end{pmatrix} = \begin{pmatrix} d_S & -\eta \\ -\mu^* & d_R \end{pmatrix}\] \end{itemize} Now we compute $V^{-1}$ and the Next Generation Matrix $K = F V^{-1}$The determinant of $V$ is: \[\det(V) = d_S d_R - \eta \mu^*\] The inverse matrix $V^{-1}$ is: \[V^{-1} = \frac{1}{d_S d_R - \eta \mu^*} \begin{pmatrix} d_R & \eta \\ \mu^* & d_S \end{pmatrix}\] Multiplying $F$ by $V^{-1}$ gives the Next Generation Matrix $K$: \[K = F V^{-1} = \frac{1}{d_S d_R - \eta \mu^*} \begin{pmatrix} \beta_S d_R & \beta_S \eta \\ \beta_R \mu^* & \beta_R d_S \end{pmatrix}\] Therefore we calculate the Reproductive Number $R_0 = \rho(K)$The reproductive number $R_0$ is the spectral radius (largest eigenvalue) of $K$.Trace of $K$: \[\text{Tr}(K) = \frac{\beta_S d_R + \beta_R d_S}{d_S d_R - \eta \mu^*}\] Determinant of $K$: \[\det(K) = \frac{\beta_S \beta_R}{d_S d_R - \eta \mu^*}\] The characteristic equation $\lambda^2 - \text{Tr}(K)\lambda + \det(K) = 0$ yields eigenvalues: \[\lambda = \frac{\text{Tr}(K) \pm \sqrt{[\text{Tr}(K)]^2 - 4\det(K)}}{2}\] Simplifying the discriminant term inside the square root: \[[\text{Tr}(K)]^2 - 4\det(K) = \frac{(\beta_S d_R + \beta_ R d_S)^2 - 4\beta_S \beta_R (d_S d_R - \eta \mu^*)}{(d_S d_R - \eta \mu^*)^2} = \frac{(\beta_S d_R - \beta_R d_S)^2 + 4 \beta_S \beta_R \eta \mu^*}{(d_S d_R - \eta \mu^*)^2}\] Taking the largest eigenvalue (with the $+$ sign), the reproductive number is: \[R_0 = \frac{(\beta_S d_R + \beta_R d_S) + \sqrt{(\beta_S d_R - \beta_R d_S)^2 + 4 \beta_S \beta_R \eta \mu^*}}{2(d_S d_R - \eta \mu^*)}\] where: \[ \beta_S = r_S + \alpha M_2^*,~~\beta_R = r_R + \alpha M_2^*,~~\mu^* = \mu_ 0 + \mu_1 M_2^* + \mu_2 C^*,~~ d_S = \delta_1 M_1^* + \kappa_1 E^* + \psi C^* + \mu^*,\] \[d_R = \delta_2 M_1^* + \kappa_2 E^* + \varepsilon \psi C^* + \eta,~ \beta_S = r_S + \alpha M _2^*,~~ \beta_R = r_R + \alpha M_2^*,~~ \mu^* = \mu_0 + \mu_1 M _2^* + \mu_2 C^*\] \[d_S = \delta_1 M_1^* + \kappa_1 E^* + \psi C^* + \mu^*,~~ \text{and}~~ d_R = \delta_2 M_1^* + \kappa_2 E^* + \varepsilon \psi C^* + \eta ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\] However, numerically the reproduction number $R_{0}=1.4467$ implying that $R_{0}>1$ using the baseline parameters, which means that cancer cell population is growing signaling a self-sustaining propagation. \subsection{Global Sensitivity Analysis}\label{app:extra3} \begin{table}[htbp] \centering \caption{Parameters with High Original S1 and Negligible Reduced S1} \label{tab:app2} \footnotesize \renewcommand{\arraystretch}{0.8} \begin{tabular}{lrrrrrrr} \hline & Original S1 & Reduced S1 & Difference S1 & Original ST & Reduced ST & Difference ST & Percentage\\ & & & & & & &Change S1 \\ Parameter & & & & & & & \\ \hline $\delta$ & $0.2057$ & $-0.0232$ & $-0.2289$ & $0.1362$ & $0.1816$ & $0.0454$ & $-111.2917$ \\ $a_2$ & $0.2026$ & $-0.0325$ & $-0.2351$ & $0.1401$ & $0.1817$ & 0.0416 & -116.0383 \\ $\mu_{M1}$ & $0.1931$ & $-0.1066$ & $-0.2997$ & $0.1540$ & $0.1766$ & $0.0226$ & $-155.2052$ \\ $k_C$ & $0.2024$ & $-0.0356$ & $-0.2380$ & $0.1577$ & $0.2210$ & $0.0633$ & $-117.5912$ \\ $K_T$ & $0.0727$ & $-0.0596$ & $-0.1324$ & $0.1954$ & $0.3712$ & $0.1758$ & $-181.9885$ \\ $\mu_E$ & $0.1534$ & $-0.0062$ & $-0.1595$ & $0.2022$ & $0.0967$ & $-0.1055$ & $-104.0237$ \\ $\mu_0$ & $0.1981$ & $-0.0203$ & $-0.2184$ & $0.2031$ & $0.1214$ & $-0.0817$ & $-110.2640$ \\ $\rho_0$ & $0.1517$ & $-0.0163$ & $-0.1680$ & $0.2044$ & $0.1005$ & $-0.1039$ & $-110.7546$ \\ $\mu_2$ & $0.1962$ & $-0.0132$ & $-0.2094$ & $0.2051$ & $0.1216$ & $-0.0836$ & $-106.7511$ \\ $\eta$ & $0.1548$ & $-0.0020$ & $-0.1568$ & $0.2060$ & $0.0941$ & $-0.1119$ & $-101.3007$ \\ $a_0$ & $0.1960$ & $-0.0126$ & $-0.2086$ & $0.2088$ & $0.1591$ & $-0.0496$ & $-106.4456$ \\ $\mu_{M2}$ & $0.2147$ & $-0.0184$ & $-0.2331$ & $0.2101$ & $0.3977$ & $0.1876$ & $-108.5705$ \\ $\delta_2$ & $0.2458$ & $-0.0313$ & $-0.2772$ & $0.2127$ & $0.1034$ & $-0.1093$ & $-112.7517$ \\ $\phi$ & $0.1926$ & $0.0046$ & $-0.1879$ & $0.2151$ & $0.1248$ & $-0.0903$ & $-97.5899$ \\ $\theta$ & $0.2166$ & $-0.0631$ & $-0.2798$ & $0.2398$ & $0.3905$ & $0.1507$ & $-129.1507$ \\ $\Lambda_E$ & $0.2046$ & $-0.0627$ & $-0.2673$ & $0.2667$ & $0.2937$ & $0.0270$ & $-130.6396$ \\ $u_0$ & $0.2096$ & $-0.0467$ & $-0.2563$ & $0.2684$ & $0.3517$ & $0.0832$ & $-122.2657$ \\ $\mu_1$ & $0.1652$ & $-0.0048$ & $-0.1699$ & $0.3014$ & $0.0939$ & $-0.2074$ & $-102.8758$ \\ $\alpha$ & $0.2636$ & $-0.1018$ & $-0.3654$ & $0.3188$ & $0.2169$ & $-0.1018$ & $-138.6018$ \\ $p_M2$ & $0.2140$ & $-0.0648$ & $-0.2789$ & $0.3308$ & $0.6401$ & $0.3093$ & $-130.2996$ \\ \hline \end{tabular} \end{table} \subsection*{Declarations} %\subsubsection*{Ethical Approval and accordance} %Not applicable. \subsubsection*{ Consent to participate} Not applicable. \subsubsection*{Consent to publish} The authors declare that they have given their consent for the publication of this manuscript, and no third-party permissions are required. \subsubsection*{Data availability} All data generated or analyzed during this study are included within the article. \subsubsection*{Competing interest} The authors declare that they have no known competing financial or personal relationships that could have appeared to influence the work reported in this study. \subsubsection*{Funding} This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors \subsubsection*{Competing interest} The authors declare that they have no known competing financial or personal relationships that could have appeared to influence the work reported in this study. %\subsubsection*{Funding} %This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors % % %\subsubsection*{Acknowledgment} %The authors are grateful to their respective institutions for the non-financial support provided during the course of this study. The views, opinions, assumptions, or any other information expressed in this article are solely those of the authors. % % % %\subsubsection*{CRediT authorship contribution statement} %\begin{itemize} %\item Conceptualization: %\item Formal analysis, Investigation, and Methodology: %\item Software: %\item Supervision: %\item Writing-original draft: %\item Writing-review and editing: %\end{itemize} % % % \subsubsection*{Declaration of generative AI and AI-assisted technologies in the writing process} Generative artificial intelligence tools were used solely for language editing and the preparation of schematic illustrations. All mathematical analyses, model development, numerical simulations, interpretation of results, and scientific conclusions were performed and verified by the authors. The authors assume full responsibility for the content of the published article. \newpage \begin{thebibliography}{99} \bibitem{ahmad2025analytical} Ahmad, W., Ullah, H., Rafiq, M., Butt, A.I.K. and Ahmad, N., 2025. Analytical and numerical investigations of optimal control techniques for managing Ebola virus disease. \textit{The European Physical Journal Plus}, 140(4), p.329. \url{https://doi.org/10.1140/epjp/s13360-025-06251-x} \bibitem{ahmad2024mathematical} Ahmad, A., Kulachi, M.O., Farman, M., Junjua, M.U.D., Bilal Riaz, M. and Riaz, S., 2024.Mathematical modeling and control of lung cancer with IL 2 cytokine and anti-PD-L1 inhibitor effects for low immune individuals. \textit{Plos one}, 19(3), p.e0299560. \url{https://doi.org/10.1371/journal.pone.0299560} \bibitem{ashrafi2022current} Ashrafi, A., Akter, Z., Modareszadeh, P., Modareszadeh, P., Berisha, E., Alemi, P.S., Chacon Castro, M.D.C., Deese, A.R. and Zhang, L., 2022. Current landscape of therapeutic resistance in lung cancer and promising strategies to overcome resistance. \textit{Cancers}, 14(19), p.4562. \url{https://doi.org/10.3390/cancers14194562} \bibitem{basak2023tumor} Basak, U., Sarkar, T., Mukherjee, S., Chakraborty, S., Dutta, A., Dutta, S., Nayak, D., Kaushik, S., Das, T. and Sa, G., 2023. Tumor-associated macrophages: an effective player of the tumor microenvironment. \textit{Frontiers in immunology}, 14, p.1295257. \url{https://doi.org/10.3389/fimmu.2023.1295257} \bibitem{bracci2014immune} Bracci, L., Schiavoni, G., Sistigu, A. and Belardelli, F., 2014. Immune-based mechanisms of cytotoxic chemotherapy: implications for the design of novel and rationale-based combined treatments against cancer. \textit{Cell Death \& Differentiation}, 21(1), pp.15-25. \url{https://doi.org/10.1038/cdd.2013.67} \bibitem{cartelle2019computational} Cartelle Gestal, M., Dedloff, M.R. and Torres-Sangiao, E., 2019. Computational health engineering applied to model infectious diseases and antimicrobial resistance spread. \textit{Applied Sciences}, 9(12), p.2486. \url{https://doi.org/10.3390/app9122486} \bibitem{chen2023macrophages} Chen, S., Saeed, A.F., Liu, Q., Jiang, Q., Xu, H., Xiao, G.G., Rao, L. and Duo, Y., 2023. Macrophages in immunoregulation and therapeutics. Signal transduction and targeted therapy. \textit{Signal Transduction and Targeted Therapy}, 8(1), p.207. \url{https://doi.org/10.1038/s41392-023-01452-1} \bibitem{donze2011robustness} Donzé, A., Fanchon, E., Gattepaille, L.M., Maler, O. and Tracqui, P., 2011. \textit{Robustness analysis and behavior discrimination in enzymatic reaction networks} PloS one, 6(9), p.e24246. \url{https://doi.org/10.1371/journal.pone.0024246} \bibitem{dunsmore2024timing} Dunsmore, G., Guo, W., Li, Z., Bejarano, D.A., Pai, R., Yang, K., Kwok, I., Tan, L., Ng, M., De La Calle Fabregat, C. and Yatim, A., 2024. Timing and location dictate monocyte fate and their transition to tumor-associated macrophages. \textbf{Science immunology}, 9(97), p.eadk3981. \url{10.1126/sciimmunol.adk3981} \bibitem{elinav2013inflammation} Elinav, E., Nowarski, R., Thaiss, C.A., Hu, B., Jin, C. and Flavell, R.A., 2013. Inflammation-induced cancer: crosstalk between tumours, immune cells and microorganisms. \textit{Nature Reviews Cancer}, 13(11), pp.759-771. \url{https://doi.org/10.1038/nrc3611} \bibitem{Eftimie2021} Eftimie, R. and Barelle, C., 2021. Mathematical investigation of innate immune responses to lung cancer: The role of macrophages with mixed phenotypes. \textit{Journal of Theoretical Biology}, 524, Article 110739. \url{https://doi.org/10.1016/j.jtbi.2021.110739} \bibitem{emran2022multidrug} Emran, T.B., Shahriar, A., Mahmud, A.R., Rahman, T., Abir, M.H., Siddiquee, M., Ahmed, H., Rahman, N., Nainu, F., Wahyudin, E. and Mitra, S., 2022. Multidrug resistance in cancer: understanding molecular mechanisms, immunoprevention and therapeutic approaches. \textit{Frontiers in Oncology}, 12, p.891652. \url{https://doi.org/10.3389/fonc.2022.891652} \bibitem{fernandez2019optimal} Fernández, L.A. and Pola, C., 2019. Optimal control problems for the Gompertz model under the Norton-Simon hypothesis in chemotherapy. \textit{Discrete \& Continuous Dynamical Systems-Series B}, 24(6), p.2577. \url{https://doi.org/10.3934/dcdsb.2018266} \bibitem{garg2024emerging} Garg, P., Malhotra, J., Kulkarni, P., Horne, D., Salgia, R. and Singhal, S.S., 2024. Emerging therapeutic strategies to overcome drug resistance in cancer cells. \textit{Cancers}, 16(13), p.2478. \url{https://doi.org/10.3390/cancers16132478} \bibitem{galete2025mathematical} Galete, M.H., Dawed, M.Y. and Obsu, L.L., 2025. Mathematical Model Dynamics of Drug‐Resistant Tuberculosis and Its Optimal Control Analysis. \textit{Journal of Applied Mathematics}, 2025(1), p.4101918. \url{https://doi.org/10.1155/jama/4101918} \bibitem{giudicianni2024variance} Giudicianni, C., Di Cicco, I., Di Nardo, A. and Greco, R., 2024. Variance-based Global Sensitivity Analysis of Surface Runoff Parameters for Hydrological Modeling of a Real Peri-urban Ungauged Basin: C. \textit{Giudicianni et al. Water Resources Management}, 38(8), pp.3007-3022. \url{https://doi.org/10.1007/s11269-024-03802-2} \bibitem{gottesman2002mechanisms} Gottesman, M.M., 2002. Mechanisms of cancer drug resistance. \textit{Annual review of medicine}, 53(1), pp.615-627. \url{https://doi.org/10.1146/annurev.med.53.082901.103929} \bibitem{hanada2014genetic} Hanada, K. and Yamaoka, Y., 2014. Genetic battle between Helicobacter pylori and humans. The mechanism underlying homologous recombination in bacteria, which can infect human cells. \textit{Microbes and infection}, 16(10), pp.833-839. \url{https://doi.org/10.1016/j.micinf.2014.08.001} \bibitem{hurvich1989regression} Hurvich, C.M. and Tsai, C.L., 1989. Regression and time series model selection in small samples. \textit{Biometrika}, 76(2), pp.297-307. \url{https://doi.org/10.1093/biomet/76.2.297} \bibitem{lassounon2026optimal} Lassounon, D., Belmiloudi, A. and Haddou, M., 2026. Optimal Control Problem with Mixed Control and State Constraints for Cancer Chemotherapy and Treatment Optimization. \textit{arXiv preprint arXiv}:2606.21765. \url{https://doi.org/10.48550/arXiv.2606.21765} \bibitem{li2025drug} Li, J., Hu, J., Yang, Y., Zhang, H., Liu, Y., Fang, Y., Qu, L., Lin, A., Luo, P., Jiang, A. and Wang, L., 2025. Drug resistance in cancer: molecular mechanisms and emerging treatment strategies. \textit{Molecular biomedicine}, 6(1), p.111. \url{https://doi.org/10.1186/s43556-025-00352-w} \bibitem{li2026advancing} Li, D., Tian, J. and Zou, Y., 2026. Advancing the understanding of the immune escape in lung cancer immunotherapy: Global trends, collaborations, and future directions. \textit{Human Vaccines \& Immunotherapeutics}, 22(1), p.2628395. \url{https://pubmed.ncbi.nlm.nih.gov/41661234/} \bibitem{li2025mitochondrial} Li, J., Zhang, L., Peng, J., Zhao, C., Li, W., Yu, Y., Huang, X., Yang, F., Deng, X., Yang, X. and Zhang, T., 2025. Mitochondrial metabolic regulation of macrophage polarization in osteomyelitis and other orthopedic disorders: mechanisms and therapeutic opportunities. Frontiers in Cell and Developmental Biology, 13, p.1604320. \bibitem{lei2023understanding} Lei, Z.N., Tian, Q., Teng, Q.X., Wurpel, J.N., Zeng, L., Pan, Y. and Chen, Z.S., 2023. Understanding and targeting resistance mechanisms in cancer. \textit{MedComm}, 4(3), p.e265. \url{https://doi.org/10.1002/mco2.265} \bibitem{li2024rabeprazole} Li, Y., Hao, J., Kong, X., Yuan, W., Shen, Y., Hui, Z. and Lu, X., 2024. Rabeprazole mitigates obesity-induced chronic inflammation and insulin resistance associated with increased M2-type macrophage polarization. \textit{Molecular Basis of Disease}, 1870(5), p.167142. \url{https://doi.org/10.1016/j.bbadis.2024.167142.} \bibitem{liu2025drug} Liu, S., Jiang, A., Tang, F., Duan, M. and Li, B., 2025. Drug-induced tolerant persisters in tumor: mechanism, vulnerability and perspective implication for clinical treatment. \textit{Molecular cancer}, 24(1), p.150. \url{https://doi.org/10.1186/s12943-025-02323-9} \bibitem{jacobs1990risk} Jacobs, I., Oram, D., Fairbanks, J., Turner, J., Frost, C. and Grudzinskas, J.G., 1990. A risk of malignancy index incorporating CA 125, ultrasound and menopausal status for the accurate preoperative diagnosis of ovarian cancer. BJOG: \textit{An International Journal of Obstetrics \& Gynaecology}, 97(10), pp.922-929. \url{https://doi.org/10.1111/j.1471-0528.1990.tb02448.x} \bibitem{kareva2015metronomic} Kareva, I., Waxman, D.J. and Klement, G.L., 2015. Metronomic chemotherapy: an attractive alternative to maximum tolerated dose therapy that can activate anti-tumor immunity and minimize therapeutic resistance. \textit{Cancer letters}, 358(2), pp.100-106. \url{https://doi.org/10.1016/j.canlet.2014.12.039} \bibitem{khalaf2021aspects} Khalaf, K., Hana, D., Chou, J.T.T., Singh, C., Mackiewicz, A. and Kaczmarek, M., 2021. Aspects of the tumor microenvironment involved in immune resistance and drug resistance. \textit{Frontiers in immunology}, 12, p.656364. \url{ https://doi.org/10.3389/fimmu.2021.656364} \bibitem{kumagai2025immunogenomic} Kumagai, S., Momoi, Y. and Nishikawa, H., 2025. Immunogenomic cancer evolution: a framework to understand cancer immunosuppression. \textit{Science Immunology}, 10(105), p.eabo5570. \url{https://doi.org/10.1126/sciimmunol.abo5570} \bibitem{khayati2023potential} Khayati, S., Dehnavi, S., Sadeghi, M., Afshari, J.T., Esmaeili, S.A. and Mohammadi, M., 2023. The potential role of miRNA in regulating macrophage polarization. Heliyon, 9(11). \textit{ Cellular and Molecular Life Sciences}, Heliyon 9(11):e21615 Published 2023 Nov. \url{https://doi.org/10.1016/j.heliyon.2023.e21615.} \bibitem{Kuznetsov1994} Kuznetsov, V.A., Makalkin, I.A., Taylor, M.A. and Perelson, A.S., 1994.Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis. \textit{Bulletin of Mathematical Biology}, 56(2), pp.295--321. \url{https://doi.org/10.1007/BF02460644} \bibitem{okosun2013impact} Okosun, K.O., Makinde, O.D. and Takaidza, I., 2013. Impact of optimal control on the treatment of HIV/AIDS and screening of unaware infectives. \textit{Applied mathematical modelling}, 37(6), pp.3802-3820. \url{https://doi.org/10.1016/j.apm.2012.08.004} \bibitem{mangalam2025ai} Mangalam, M., 2025. AI-driven dynamic grouping for adaptive clinical trials: Rethinking randomization in precision medicine. \textit{Artificial Intelligence in Medicine}, p.103272. \url{10.1016/j.artmed.2025.103272} \bibitem{Mantovani2008} Mantovani, A., Allavena, P., Sica, A. and Balkwill, F., 2008. Cancer-related inflammation. \textit{Nature}, 454(7203), pp.436--444. \url{https://doi.org/10.1038/nature07205} \bibitem{ma2024comprehensive} Ma, C. and Gurkan-Cavusoglu, E., 2024. A comprehensive review of computational cell cycle models in guiding cancer treatment strategies. \textit{NPJ Systems Biology and Applications}, 10(1), p.71. \url{https://doi.org/10.1038/s41540-024-00397-7} \bibitem{mohammad2021mathematical} Mohammad Mirzaei, N., Su, S., Sofia, D., Hegarty, M., Abdel-Rahman, M.H., Asadpoure, A., Cebulla, C.M., Chang, Y.H., Hao, W., Jackson, P.R. and Lee, A.V., 2021. A mathematical model of breast tumor progression based on immune infiltration. \textit{Journal of Personalized Medicine}, 11(10), p.1031. \url{https://doi.org/10.3390/jpm11101031} \bibitem{Murray2011} Murray, P.J. and Wynn, T.A., 2011. Protective and pathogenic functions of macrophage subsets. \textit{Nature Reviews Immunology}, 11(11), pp.723--737. \url{https://doi.org/10.1038/nri3073} \bibitem{nam2024harnessing} Nam, Y., Kim, J., Jung, S.H., Woerner, J., Suh, E.H., Lee, D.G., Shivakumar, M., Lee, M.E. and Kim, D., 2024. Harnessing artificial intelligence in multimodal omics data integration: paving the path for the next frontier in precision medicine. \textit{Annual review of biomedical data science}, 7(1), pp.225-250. \url{https://doi.org/10.1146/annurev-biodatasci-102523-103801} \bibitem{neath2012bayesian} Neath, A.A. and Cavanaugh, J.E., 2012. The Bayesian information criterion: background, derivation, and applications. Wiley Interdisciplinary Reviews: \textit{Computational Statistics}, 4(2), pp.199-203. \url{https://doi.org/10.1002/wics.199} \bibitem{Noy2014} Noy, R. and Pollard, J.W., 2014. Tumor-associated macrophages: From mechanisms to therapy. \textit{Immunity}, 41(1), pp.49--61. \url{https://doi.org/10.1016/j.immuni.2014.06.010} \bibitem{pacheco2017cancer} Pacheco, E.S.M., 2017. Cancer therapies based on Optimal Control methods. \bibitem{pradhan1995steady} Pradhan, A., Pal, P., Durocher, G., Villeneuve, L., Balassy, A., Babai, F., Gaboury, L. and Blanchard, L., 1995. Steady state and time-resolved fluorescence properties of metastatic and non-metastatic malignant cells from different species. \textit{Journal of Photochemistry and Photobiology B: Biology}, 31(3), pp.101-112. \url{https://doi.org/10.1016/1011-1344(95)07178-4} \bibitem{ramos2021battling} Ramos, A., Sadeghi, S. and Tabatabaeian, H., 2021. Battling chemoresistance in cancer: root causes and strategies to uproot them. \textit{International journal of molecular sciences}, 22(17), p.9451. \url{https://www.google.com/search?q=https://doi.org/10.3930/ijms22179451} \bibitem{siannis2005sensitivity} Siannis, F., Copas, J. and Lu, G., 2005. Sensitivity analysis for informative censoring in parametric survival models. \textit{Biostatistics}, 6(1), pp.77-91. \url{https://doi.org/10.1093/biostatistics/kxh019} \bibitem{solinas2009tumor} Solinas, G., Germano, G., Mantovani, A. and Allavena, P., 2009. Tumor-associated macrophages (TAM) as major players of the cancer-related inflammation. \textit{Journal of Leukocyte Biology}, Volume 86, Issue 5, Nov 2009 Pages 1065–1073. \url{https://doi.org/10.1189/jlb.0609385} \bibitem{touray2025combating} Touray, B.J., 2025. Combating Persistent and Emerging Respiratory Pathogens Through Innovative Nanovaccines. \textit{The University of Wisconsin-Madison}. \url{https://www.proquest.com/openview/729fe7ffcb8c96177af84d908bd70136/1?pq-origsite=gscholar&cbl=18750&diss=y} \bibitem{tufail2025immune} Tufail, M., Jiang, C.H. and Li, N., 2025. Immune evasion in cancer: mechanisms and cutting-edge therapeutic approaches. \textit{Signal transduction and targeted therapy}, 10(1), p.227. \url{https://doi.org/10.1038/s41392-025-02280-1} \bibitem{Enderling2014} Enderling, H. and Chaplain, M.A.J., 2014. Mathematical modeling of tumor growth and treatment. \textit{Current Pharmaceutical Design}, 20(30), pp.4934--4940. \url{https://doi.org/10.2174/1381612819666131125150434} \bibitem{Foo2012} Foo, J. and Michor, F., 2012. Evolution of resistance to targeted anti-cancer therapies during continuous and pulsed administration strategies. \textit{PLoS Computational Biology}, 8(11), e1002557. \url{https://doi.org/10.1371/journal.pcbi.1002557} \bibitem{dePillis2005} de Pillis, L.G., Radunskaya, A.E. and Wiseman, C.L., 2005. A validated mathematical model of cell-mediated immune response to tumor growth. \textit{Cancer Research}, 65(17), pp.7950--7958. \url{https://doi.org/10.1158/0008-5472.CAN-05-0564} \bibitem{jamal2017tracking} Jamal-Hanjani, M., Wilson, G.A., McGranahan, N., Birkbak, N.J., Watkins, T.B., Veeriah, S., Shafi, S., Johnson, D.H., Mitter, R., Rosenthal, R. and Salm, M., 2017. Tracking the evolution of non–small-cell lung cancer. \textit{New England Journal of Medicine}, 376(22), pp.2109-2121. \url{https://www.nejm.org/doi/full/10.1056/NEJMoa1616288} \bibitem{jeyachandran2025invertebrate} Jeyachandran, S. and Maran, B.A.V., 2025. Invertebrate Immunology. \textit{Springer Nature Singapore}. \url{https://link.springer.com/book/10.1007/978-981-95-1549-3} \bibitem{he2026macrophage} He, W., Xu, J. and Li, X., 2026. Macrophage polarization in inflammatory regulation: Molecular mechanisms, therapeutic targets, and translational challenges. Cellular and Molecular Life Sciences, 83(1), p.164. \textit{Cellular and Molecular Life Sciences} \url{https://link.springer.com/article/10.1007/s00018-025-06041-9#citeas} \bibitem{injerd2017mathematical} Injerd, R. and Turian, E., 2017. Mathematical Modeling of Non-Small Cell Lung Cancer Response to Therapy (No. 17-0922). \textit{Technical Report}. \url{https://www.neiu.edu/sites/default/files/inline_files/Russell_TechnicalReport_Apr2018.pdf} \bibitem{Gatenby2009} Gatenby, R.A., Silva, A.S., Gillies, R.J. and Frieden, B.R., 2009. Adaptive therapy. \textit{Cancer Research}, 69(11), pp.4894--4903. \url{https://doi.org/10.1158/0008-5472.CAN-08-3658} \bibitem{Coldman1983} Coldman, A.J. and Goldie, J.H., 1983. A stochastic model for the origin and treatment of tumors containing drug-resistant cells. \textit{Bulletin of Mathematical Biology}, 45(5), pp.739--757. \url{https://doi.org/10.1016/S0092-8240(83)80012-7} \bibitem{Coddington1955} E. A. Coddington and N. Levinson,\textit{Theory of Ordinary Differential Equations},McGraw--Hill, New York, 1955. \bibitem{miller2018stable} Miller, J. and Hardt, M., 2018. \textit{Stable recurrent models} arXiv preprint arXiv:1805.10369. \url{https://arxiv.org/abs/1805.10369} \bibitem{oliveira2024overview} Oliveira, M., Antunes, W., Mota, S., Madureira-Carvalho, Á., Dinis-Oliveira, R.J. and Dias da Silva, D., 2024. An overview of the recent advances in antimicrobial resistance. \textit{Microorganisms}, 12(9), p.1920. \url{ https://doi.org/10.3390/microorganisms12091920} \bibitem{Perko2001} L. Perko,\textit{Differential Equations and Dynamical Systems},3rd ed., Springer, New York, 2001. \bibitem{van2002reproduction} Van den Driessche, P. and Watmough, J., 2002. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. \textit{Mathematical biosciences}, 180(1-2), pp.29-48. \url{https://doi.org/10.1016/S0025-5564(02)00108-6} \bibitem{van2018molecular} Van Dalen, F.J., Van Stevendaal, M.H., Fennemann, F.L., Verdoes, M. and Ilina, O., 2018. \textit{Molecular repolarisation of tumour-associated macrophages. Molecules}, 24(1), p.9. \url{https://doi.org/10.3390/molecules24010009} \bibitem{wang2025lung} Wang, Z., Guo, H., Song, Y., Wang, A., Yan, Y., Ma, L. and Liu, B., 2025. Lung cancer tumor immune microenvironment: analyzing immune escape mechanisms and exploring emerging therapeutic targets. \textit{Frontiers in immunology}, 16, p.1597686. \url{https://doi.org/10.3389/fimmu.2025.1597686} \bibitem{xiong2025mathematical} Xiong, Z., Xia, Y., Xue, L. and Lei, J., 2025. Mathematical Modelling and Optimization of Medication Regimens for Combination Immunotherapy of Breast Cancer: Z. Xiong et al. \textit{Bulletin of Mathematical Biology}, 87(7), p.88. \url{https://doi.org/10.1007/s11538-025-01459-5} \bibitem{yan2025progress} Yan, W., Li, X., Xu, S., Pan, H., Wang, B., Zeng, Z. and Shi, Y., 2025. Progress in research on macrophage polarization mechanisms and targeted therapies in Staphylococcus aureus infections. \textit{Infection and Drug Resistance}, pp.6655-6671. \url{https://doi.org/10.2147/IDR.S551067} \bibitem{young2024data} Young, A.L., Oxtoby, N.P., Garbarino, S., Fox, N.C., Barkhof, F., Schott, J.M. and Alexander, D.C., 2024. Data-driven modelling of neurodegenerative disease progression: thinking outside the black box. \textit{Nature Reviews Neuroscience}, 25(2), pp.111-130. \url{https://doi.org/10.1038/s41583-023-00779-6} \bibitem{zaider2011tumor} Zaider, M. and Hanin, L., 2011. Tumor control probability in radiation treatment. \textit{Medical physics}, 38(2), pp.574-583. \url{https://doi.org/10.1118/1.3521406} \bibitem{zhang2025mathematical} Zhang, H. and Li, C., 2025. Mathematical modelling of tumor-immune interactions in breast cancer. \textit{Journal of Theoretical Biology}, p.112310. \url{https://doi.org/10.1016/j.jtbi.2025.112310} \bibitem{zhou2026bacteriophage} Zhou, Z., Fu, H., Li, M., Han, Z., Wu, Z., Fan, H., Shen, N. and Zheng, J., 2026. \textit{Bacteriophage therapy: current strategies and future perspectives}. MedComm, 7(3), p.e70645. \url{ https://doi.org/10.1002/mco2.70645} \end{thebibliography} \end{document}