Cancer progression is frequently driven by mutations in the tumor suppressor gene TP53. These mutations cause the protein to lose its normal growth-inhibitory functions and acquire new properties that promote tumor development and immune evasion. Recent experimental studies show that mutant p53 (mutp53) suppresses the immune system by binding to TANK-binding kinase 1 (TBK1). This interaction blocks the formation of the Stimulator of Interferon Genes (STING)-TBK1-interferon regulatory factor 3 (IRF3) signaling complex, which is essential for inducing Type I interferon (IFN-I) production. Consequently, the anti-tumor immune response is weakened. Despite increasing experimental evidence for the mutp53-TBK1 interaction, the quantitative impact of this axis on tumor-immune dynamics remains poorly understood. In particular, it is unclear how modulation of mutp53 or restoration of TBK1 activity influences long-term tumor control at the systems level. To investigate the potential benefits of mutp53-TBK1-based therapeutic interventions, we developed a mathematical model based on a system of ordinary differential equations (ODEs). This proposed model captures the dynamic interplay between tumor growth, immune effector cells, IFN-I, and the intracellular STING-TBK1-IRF3 signaling axis regulated by mutp53. Parameter uncertainty was explored using Latin Hypercube Sampling (LHS), and the resulting model outputs were analyzed using Partial Rank Correlation Coefficients (PRCC) to evaluate parameter significance and model robustness. Numerical simulations yield key predictive insights into tumor-immune dynamics. First, in the absence of treatment, the model indicates that mutp53-mediated suppression of innate signaling drives immune escape and sustained tumor progression. Second, under simulated therapeutic conditions, the results suggest a profound divergence in treatment robustness: shRNA-mediated knockdown of mutp53 is predicted to yield highly consistent tumor suppression regardless of the baseline immune state. Conversely, the efficacy of ectopic TBK1 overexpression relies heavily on existing immune recruitment, with its impact diminishing significantly in immunologically 'cold' environments. Finally, the model predicts that combining both interventions produces a strong synergistic antitumor effect, achieving the greatest reduction in tumor burden. These theoretical results, which quantitatively align with recent experimental observations, underscore the critical need to simultaneously target upstream mutp53 oncogenic signaling and augment downstream TBK1 induction to successfully re-establish effective immune surveillance.
This work examines a nonlinear stochastic SIRS epidemic model evolving in a randomly changing environment described by a finite-state Markov chain. The transmission mechanism incorporates regime-dependent nonlinear incidence rates, where the contact interaction between susceptible and infectious individuals is modeled by the term [Formula: see text] . This switching nonlinearity allows the model to capture varying environmental effects and heterogeneous transmission patterns more accurately, thereby providing a more realistic description of epidemic dynamics. To the best of our knowledge, the sufficient criteria governing the persistence and extinction of stochastic SIRS models with transmission rate exponents governed by Markovian switching have not yet been established in the existing literature. The principal contribution of the present study is the derivation of the rigorous sufficient conditions that characterize both the extinction and the long-term persistence of disease dynamics. Specifically, a threshold parameter Λ, expressed in terms of the switching exponents ρξ(t) and ζξ(t), is derived. That is, if Λ > 0, the disease exhibits strong stochastic persistence; conversely, if Λ < 0, the disease-free equilibrium state becomes globally asymptotically stable in probability, leading to eventual disease extinction. In the special case where there is no regime switching and ρξ(t)=ζξ(t)=1, our model recovers the classical threshold found in the literature. To support and validate the theoretical findings, numerical simulations are provided to demonstrate the dynamical behavior of the model under different environmental regimes.
The Wilson-Cowan neural mass model groups excitatory and inhibitory neurons and models their communication through a system of ODEs. The neural firing rate function is popularly approximated by a smooth sigmoid function and a discontinuous Heaviside function. The nondifferentiability of the latter function invalidates standard methods that require derivative information, such as stability and sensitivity theory. In this article, we directly analyze a Wilson-Cowan model with a piecewise linear firing rate function (modeling the "all-or-nothing" action potential) using a relatively new tool from generalized derivatives theory called the lexicographic directional derivative. Our contributions include establishing well-posedness of the nonsmooth Wilson-Cowan model and analyzing its stability and parametric sensitivities. In the smooth Wilson-Cowan model the limit cycle, which corresponds to spiking behavior, is globally attractive in the physically-meaningful domain, while the limit cycle is only locally attractive in the nonsmooth model because a locally attractive trivial equilibrium also exists in that case. A local and nonlocal sensitivity analysis of the nonsmooth model uncovers that the inhibitory neuron time scale is the most influential parameter in the nonsmooth case, dominating the effects from the other parameters. This differs significantly from the smooth model, which we also observe to be overall less sensitive to its parameters than the nonsmooth model, despite the state variable solutions in each model following similar spiking behaviors.
Cancer progression emerges from a dynamical interplay between neoplastic growth and immune surveillance, yet most mathematical models reduce this interaction to logistic kinetics or single-threshold Allee effects, overlooking the layered barrier structure that governs tumor establishment. We introduce a three-dimensional immuno-oncology model coupling immune effector cells, cancer cells, and immunotherapy, which we reduce via singular perturbation and scale to a dimensionless planar system. The key novelty lies in endowing both compartments with generalized Allee functions: a single threshold for immune recruitment and a dual hyper-Allee law for tumor growth, reflecting cooperative, autocrine, and threshold-dependent processes at low cell densities. Through bifurcation analysis and numerical continuation in the clinically motivated parameter plane of immune loss versus cytotoxic efficacy, we uncover a global atlas of qualitatively distinct regimes. The organizing center is a cusp-type degenerate Bogdanov-Takens bifurcation of codimension three, from which saddle-node, Hopf, homoclinic, and limit-point-of-cycle curves emanate. Our results provide a quantitative scaffold for the three Es of cancer immunoediting: elimination, equilibrium, and escape, and suggest that patient stratification and therapeutic timing should be viewed as navigation problems in a high-dimensional threshold landscape.
This work investigates the dynamics of ionic flows through a classical Poisson-Nernst-Planck model incorporating two oppositely charged ion species and small permanent charges. By integrating boundary layer effects into current-voltage relations, we analyze how ionic flow dynamics are influenced by the interplay of physical parameters including permanent charge, channel geometry, and ion diffusion coefficients. Within a simplified modeling framework, we identify several critical potentials that offer mechanistic insights into electrodiffusion processes in biological channels. Extending the work presented in [Membranes 2023, 13, 131], we examine the sequence of these critical potentials and characterize how small permanent charges affect ion flux within the boundary layer. Numerical simulations are conducted to illustrate the analytical results intuitively.
To understand how extracellular potassium ions influence both anti-tumor immunity and virotherapy in cancer progression, we construct two mathematical models informed by our experimental results. For the tumor immune system, our analysis shows a high concentration of extracellular potassium ions diminishes the killing rate of tumor cells by immune cells. Our model confirms that the stimulation coefficient of immune cells by tumor cells remains a crucial parameter in the presence of extracellular potassium ions, which largely controls the overall tumor growth. For virotherapy, we obtain a formula for the basic viral reproduction number which combines several parameters including potassium ions. When this number is greater than one, virotherapy achieves some partial success, where the tumor load is an increasing function of the potassium ion concentration. Therefore, the tumor load is reduced if the potassium ion concentration is lowered. The delay parameter of the viral lytic cycle also affects the basic reproduction number, and we find a critical delay time which determines when the reproduction number is greater than one. We performed some numerical analysis. A two-parameter bifurcation analysis reveals a positive correlation between viral burst size and the potassium-ion absorption rate, suggesting that higher absorption rates can enhance the success of virotherapy.
This paper develops a rigorous mathematical framework capable of linking invasion thresholds, trait-dependent fitness, and evolutionary stability. We propose an integro-differential within-host HIV-1 model in which the level of antiretroviral resistance is represented as a continuous (quantitative) trait structuring infected cells and viral particles. The model incorporates key biological mechanisms, including intracellular delay, cytotoxic T lymphocyte (CTL) immune responses, and both virus-to-cell and cell-to-cell transmission pathways. From a theoretical perspective, we establish the well-posedness of the model and prove its dissipativity and asymptotic compactness. Using perturbation and spectral methods, we characterize the basic reproduction number, R0, as the spectral radius of an associated next-generation operator. We further derive a direct connection between R0 and a resistance-dependent fitness function Θ, thereby linking epidemiological invasion criteria to the adaptive landscape governing resistance evolution. We prove that the infection-free equilibrium is globally asymptotically stable when R0<1, whereas uniform persistence occurs when R0>1. Our analysis reveals a fundamental evolutionary principle: viral persistence is a necessary condition for evolutionary selection. In particular, the sign of maxxΘ(x)-1 determines whether adaptive evolution can occur, while the shape of Θ determines the location of evolutionary attractors. Numerical simulations further highlight three qualitatively distinct regimes. When treatment suppresses the maximal invasion fitness below unity, viral extinction occurs before adaptive structuring can emerge. When treatment is only partially effective, viral persistence coexists with directional selection toward resistant phenotypes, generating stable evolutionary attractors at elevated resistance levels. In contrast, in the absence of therapy, resistance-associated fitness costs dominate and selection favors low-resistance strains. Finally, we show that the long-term evolutionary outcome depends not only on the fitness landscape but also on the structure of the mutation process. Under symmetric mutation kernels, evolutionary attractors coincide with fitness optima, whereas directional mutation biases can substantially shift the dominant phenotype away from the fitness-maximizing resistance level. These results demonstrate how therapeutic pressure and mutation jointly reshape the adaptive landscape and determine the emergence, persistence, and evolutionary endpoint of drug-resistant HIV populations.
During prolonged epidemics, public adherence to protective measures often wanes due to behavioral fatigue-a dynamic psychological process typically overlooked in existing models. We propose a coupled Unaware-Positive-Negative-Unaware and Susceptible-Infected-Recovered-Susceptible (UPNU-SIRS) contagion model on a two-layer multiplex network, where both layers are constructed as 2-simplicial complexes to capture higher-order group interactions. The core innovation is an emotion-threshold mechanism: individuals accumulate negative emotions, and exceeding a heterogeneous personal threshold dynamically increases the probability of behavioral regression (e.g., from active to passive protection) and reduces immunity. Using the Microscopic Markov Chain Approach and extensive Monte Carlo simulations, we find that this mechanism fundamentally drives secondary outbreaks and sustains elevated endemic levels. Critically, while the emotion-threshold mechanism alone cannot induce bistability, its synergy with higher-order interactions in the disease layer dictates the extent of the bistable regime and the outbreak threshold-offering a leverage point where enhancing emotional resilience, reducing immunity-weakening effects, and ensuring timely, trusted media communication can collectively suppress both initial and secondary epidemic peaks. By integrating socio-psychological feedback with higher-order network structures, our framework provides new insights into epidemic persistence and targeted interventions against pandemic fatigue.
The Philadelphia chromosome-negative myeloproliferative neoplasms (MPNs) are a group of haematological malignancies triggered by a driver mutation, most commonly the JAK2 V617F mutation, which is acquired in a haematopoietic stem cell. The diseases are characterised by an overproduction of myeloid cells and may result in severe complications such as thrombosis, myelofibrotic transition, and potentially progression to acute myeloid leukaemia. Treatment with interferon-α (IFN-α) can potentially deplete the disease-driving malignant stem cell population, thereby leading to long-term remission. We extend a mechanistic compartmental differential equation model of MPN progression to account for the effects of IFN-α treatment. In agreement with experimental mouse studies, we assume that IFN-α acts through effects on malignant stem cell differentiation and malignant progenitor and precursor cell apoptosis. We use a hierarchical Bayesian inference procedure to infer the drug response parameters of the model based on measurements of the JAK2 V617F variant allele frequency (VAF) for N=56 patients from the Danish DALIAH study, both on the individual and the population level, with 14 patients left out for subsequent testing. The model estimates are found to agree with data. The drug response parameters are found to act synergistically on the reduction of JAK2 VAF, with the model being able to capture qualitatively different types of treatment responses. Using the inferred information about the population distribution of drug response parameters is found to improve predictions compared to predictions without prior knowledge on a testing cohort. Such predictions may aid clinical decision-making regarding IFN-α treatment in patients with MPNs.
Notch-Delta signaling is a juxtacrine cell-fate pathway that drives lateral inhibition, whereby neighboring cells adopt alternating Sender (high Delta/low Notch) and Receiver (high Notch/low Delta) states. This bistability is captured by the single-cell Notch-Delta circuit model with external Notch and Delta inputs introduced by Boareto et al. [1], which underlies many computational frameworks. The model incorporates negative feedback in Delta production and cis-inhibition of Notch by Delta. In this work, a rigorous mathematical analysis of the model is presented. It is shown that, generically in parameter space, the system admits at most three equilibria. For Hill coefficients equal to one and two, sufficient conditions for the existence of three equilibria are derived. When exactly three distinct nondegenerate equilibria exist and the Hill coefficient lies in 0<p<3+8, bistability is established. The quasi-steady-state (QSS) reduction obtained by eliminating the Notch intracellular domain (NICD) exhibits identical bistable behavior to the full system over the same Hill coefficient range, independent of the time-scale separation parameter. When the time-scale separation parameter approaches zero, the slow variables (Notch and Delta) of the full model converge uniformly in time to those of the QSS model, whereas the fast variable (NICD) converges after a brief initial layer. Bifurcation analyses, together with a randomized parameter sweep,reveal that the monostable Receiver and Sender phenotypes are robust and dominant, whereas the bistable dual phenotypes (Receiver/Sender, Receiver/Receiver, Sender/Sender) are rare and sensitive to parameter changes when the parameter ranges are broad. Finally, a stability analysis framework based on Determinant Reduction Theorem 1 is introduced, substantially simplifying the analysis and accommodating real-valued Hill coefficients.
Malignant melanoma is an aggressive skin cancer with limited responsiveness to traditional therapies. Notably, the combination of dendritic cell (DC) vaccines with anti-programmed cell death protein 1 (anti-PD-1) therapy has shown stronger clinical potential than conventional approaches. Understanding the tumor-immune interplay is essential for optimizing melanoma immunotherapy strategies. In this paper, we formulate a melanoma-specific tumor-immune interaction model of tumor cells (TCs), DCs, and effector CD8+ T cells (ECs). A key threshold value is identified to characterize tumor growth. Using this threshold, we determine the conditions for tumor-free and tumorous equilibria, consistent with cancer immunoediting theory. Furthermore, bifurcation analysis indicates that the model exhibits oscillatory behavior under certain conditions. Sensitivity and parameter heterogeneity analyses reveal that tumor burden is mainly regulated by the intrinsic tumor growth and immune activation rate. Moreover, to better reflect physiological realism, we extend the model to a time-delayed system by incorporating a constant delay for DC-to-EC activation. Analytical and numerical results demonstrate a supercritical Hopf bifurcation at a critical delay τ0 ≈ 4.68 days, leading to stable periodic solutions. Finally, an optimal control framework is proposed to design DC vaccines and anti-PD-1 injection protocols. Compared with the constant dosing strategy, optimal control achieves enhanced tumor suppression for the same total treatment intensity. This work elucidates the dynamical mechanisms of melanoma-immune interactions and establishes a theoretical foundation for personalized combination immunotherapies.
We present a Laplace transform approach for explicitly solving linear delay differential systems with multiple discrete delays by applying the Cauchy residue theorem. This method enables direct determination of the stability of the trivial solution when delays are relatively small. Its efficacy is illustrated through two nonlinear models with two and three delays, respectively, for which explicit solutions and stability criteria are obtained. The approach offers two key advantages: (i) analytic solutions are obtained with less effort than the method of steps, and (ii) a small hyper-tetrahedron region in the delay parameter space can be identified in which the trivial solution is asymptotically stable. Furthermore, the results can be combined with existing theory, such as Lemma 3.10 in [1], to establish conditions for Hopf bifurcation.(Dedicated to Professor Shigui Ruan on the occasion of his 60th birthday)
Atrial fibrillation (AF) is the most prevalent sustained cardiac arrhythmia and is strongly associated with electrical and structural remodeling. Fibrosis, a hallmark of structural remodeling, alters myocardial conduction through architectural disruption and electrotonic interactions between cardiomyocytes and fibroblasts. However, the combined effects of fibrotic texture, density, and cellular coupling on reentrant dynamics remain incompletely understood. In this work, a computational model integrating electrical and structural remodeling is used to investigate how different fibrotic architectures and degrees of heterogeneity influence reentrant mechanisms during AF. A two-dimensional atrial tissue model was implemented using a complex-order monodomain formulation to represent structural heterogeneities. Electrotonic coupling between cardiomyocytes and fibroblasts was incorporated using biophysically detailed ionic models. Compact, diffuse, and patchy fibrosis textures were simulated, including scenarios representing clinical Utah classification fibrosis densities. Reentrant activity was induced using an S1-S2 stimulation protocol and analyzed through phase singularity tracking and virtual electrograms. Simulations showed that reentrant waves anchor to fibrotic regions, with dynamics dependent on fibrotic texture, density, and complex derivative order. Increased heterogeneity and fibroblast coupling modulated dominant frequency, wave stability, and phase singularity trajectories. Quantitative analysis further showed that the imaginary part of the complex order increased the variability of activation dynamics while preserving the spatial organization of reentrant trajectories within each fibrosis configuration. Diffuse and patchy fibrosis produced heterogeneous reentry patterns, while compact fibrosis favored stable macroreentries under high fibroblast coupling.
Externally induced overcompensation can generate paradoxical population responses, including hydra and hormetic effects, in which external disturbances increase equilibrium population size or produce nonmonotone response curves. In this study, we formulate a single-population discrete-time model with heterogeneous harvesting and alternative placements of compensatory feedback. The model divides the population update into two parts: individuals or biomass that persist directly to the next generation, and the remaining component that contributes through density-dependent recruitment. This formulation allows us to compare how different compensatory assumptions affect equilibrium responses. When the positive equilibrium is locally stable, we derive threshold conditions for hydra and hormetic effects and show that their occurrence depends on where compensation is introduced in the population update. We also examine the local stability and bifurcation behavior of the positive equilibrium as harvesting intensity changes. The equilibrium formulas are further fitted to two published nonmonotone response datasets. The fitted curves show good agreement with the observed U-shaped and inverted U-shaped patterns.
Predators rarely forage with fixed efficiency; they cooperate, adjust their search effort to prey availability, and impose non-consumptive effects such as fear, while prey often respond through mutualistic associations with non-prey partners. Yet how these adaptive feedbacks collectively shape ecosystem resilience, critical transitions, and collapse in mutualistic communities remains unresolved within a unified theoretical framework. To address this gap, we develop a community-level mathematical model integrating prey-non-prey mutualism, cooperative hunting, prey-dependent predator search efficiency, and fear-mediated reproductive suppression, and investigate how these interacting mechanisms influence ecological stability under environmental variability. Our results reveal that cooperative hunting fundamentally reorganizes the ecosystem’s tipping structure. Increasing cooperation drives catastrophic regime shifts through saddle-node bifurcations, promotes oscillatory coexistence, and generates homoclinic transitions absent from earlier mutualistic predator-prey models. Mutualistic support initially enhances resilience and species persistence but, beyond a critical threshold, induces tri-stable dynamics consistent with the paradox of enrichment. Contrary to the conventional view of fear as solely detrimental, moderate fear stabilizes coexistence, whereas its interaction with strong predator cooperation triggers abrupt transitions between alternative ecological states. Furthermore, two-parameter bifurcation analysis uncovers a remarkably rich dynamical landscape containing Bogdanov-Takens, generalized Hopf, cusp, and homoclinic structures, exposing multiple pathways to ecosystem collapse and recovery that remained hidden in previous studies lacking adaptive predator behavior. Under environmental stochasticity, noise-induced switching emerges well before deterministic tipping thresholds are reached, while increasing standard deviation and lag-1 autocorrelation provide reliable early-warning signals of impending collapse. Together, these findings identify adaptive predator behavior as a key determinant of resilience, tipping dynamics, and collapse predictability in mutualistic ecosystems facing environmental change.
We study how a social insect colony reallocates workers between local entrance defence and other colony tasks during a short disturbance. To make this trade-off explicit, we extend a patrol-recruit framework by separating workers in a waiting pool from workers engaged in non-defence tasks such as foraging and nest work. On a behavioural time scale short relative to colony demography, we formulate a four-class model for patrollers, alarmed recruiters, waiting workers, and other-task workers, and reduce it, via workforce conservation, to a three-dimensional system. We show that all biologically feasible solutions remain bounded and that the colony cannot lose its waiting pool or entrance patrols through internal dynamics alone. The reduced system always admits a non-recruiting equilibrium and can admit at most two recruiting equilibria. Their existence depends on colony size and a scalar balance relation for the waiting pool, while recruiter invasion occurs only above a colony-size threshold. We further show that at most one recruiting equilibrium can be stable, and that the candidate stable branch is selected by a simple comparison between recruitment efficiency and turnover between waiting and other-task workers. These results identify parameter ranges with local bistability between the non-recruiting equilibrium and a recruiting equilibrium and distinguish lean from buffered defence configurations. Numerical simulations agree with the analysis and indicate convergence to steady worker-allocation patterns rather than sustained oscillations.
Traditionally employed as an ecological tool to prevent wildfires, controlled burns may also reduce the populations of vectors like ticks, thereby helping to manage vector-borne diseases. In this study, we analyze an impulsive ordinary differential equation (ODE) model that incorporates fire-driven tick-host population dynamics and the transmission of a tick-borne disease with co-feeding pathway. We assess the effects of fire frequency and intensity on tick dynamics and the transmission of tick-borne diseases. Additionally, we examine the impact of fire control when disease transmission occurs not only through systemic infection but also via co-feeding. Our findings highlight that fire control outcomes depend not only on the efficacy and periodicity of interventions, but also on the host-vector population structure, as well as the epidemiological characteristics of the disease, such as its transmission mode.
We analyze discrete-time Leslie-Gower competition models featuring Beverton-Holt survivability and interspecific mating to explore interactions between two Aedes mosquito species. By incorporating larval resource competition and asymmetric reproductive interference (satyrization), we identify conditions for equilibrium existence and local stability. We specifically examine how the mating competitiveness of Ae. albopictus and its interference with Ae. aegyptidrive competitive displacement or coexistence of the two mosquito species. Our results demonstrate that satyrization significantly alters classical competition outcomes, leading to displacement or coexistence depending on parameter regimes and initial population sizes. These findings, supported by numerical simulations, provide a robust theoretical framework for understanding how reproductive interference shapes mosquito population dynamics.
The effect of dispersal in a fragmented habitat is often studied using temporally homogeneous metapopulation models. However, natural environments are characterized by temporal heterogeneity, where factors like resource availability and temperature fluctuate over time. This study extends the classic discrete-time two-patch model by incorporating temporal heterogeneity in model parameters to investigate the effect of environmental fluctuations on total biomass. Our analysis reveals that while fluctuations preserve some established principles, such as the always-beneficial effect of connecting two source patches with identical intraspecific competition, they also generate novel and unexpected dynamics. We show that fluctuations generate new response patterns of total biomass to an increase of dispersal. Furthermore, we demonstrate that fluctuations can alter the qualitative response of connectivity, creating situations where connecting two patches, each of which would suffer from connectivity in a static environment, can become beneficial under periodic conditions. These results underscore the critical importance of incorporating temporal heterogeneity to accurately predict the consequences of habitat fragmentation and to design effective conservation strategies.
