Heterogeneous immune recovery after viral response through a dynamical model of feedback-driven persistence and clearance

Abstract

Viral infections trigger complex immune responses with heterogeneous outcomes shaped by nonlinear feedback. A system of ordinary differential equations is developed to investigate immune response dynamics during viral infection, incorporating six modules, including viral load, innate immunity, cellular immunity, humoral immunity, immune suppression, and IL-6. Bifurcation analysis reveals that under continuous viral exposure, when viral clearance rate and intrinsic viral death rate satisfy specific conditions, the system exhibits up to five stable equilibria. This indicates that different health and disease states may coexist depending on initial conditions, while severe inflammation mainly arises from strong activation of cellular immunity, highlighting the complexity of immune responses. Simulations of finite-time viral exposure demonstrate multi-timescale recovery characteristics, with viral load and IL-6 levels declining rapidly, while humoral immune activation and immunosuppression show delayed and sustained patterns. Furthermore, analysis of infectious period and disease duration also indicates that during transition from early acute response to chronic disease, viral replication rate plays a critical role, while immune response intensity is sensitive to both viral clearance and immune self-activation. Subsystem analysis identifies the three-component subsystem of viral load, innate immunity, and cellular immunity as core drivers of bistability and oscillations, while humoral immunity, immune suppression, and IL-6 primarily modulate response amplitude and timing. This work establishes a theoretical framework for analyzing immune response and chronic risks through feedback dynamical modelling, providing insights for intervention strategies.

1 Introduction

The host immune response triggered by viral infection is a highly complex and dynamic process, involving nonlinear and coupled interactions among multiple cytokines and regulatory mechanisms (1–3). These immune interactions determine the course and severity of diseases and profoundly influence immune homeostasis and tissue repair during recovery (4), which ultimately significantly shape the dynamics of epidemic transmission (5–7). Classical viral dynamics models have successfully captured the basic relationship between viral load and cell infection (8, 9), while an integrative characterization is required for the overall regulatory effects of immunity on viral dynamics (8) given the high complexity of feedback regulation within the immune system. The immune system is increasingly regarded as a nonlinear complex adaptive network shaped by intertwined activation and suppression pathways with multistability, thresholds, and oscillatory responses (10). Recently, increasing evidence has shown that multistability, sustained oscillations, and sensitivity to initial conditions are key dynamical features underlying the heterogeneity of post-infection immune states (11). For example, a COVID-19 immunity-based model demonstrated that variations in the timing and strength of innate, cellular, and inflammatory responses can generate markedly divergent clinical trajectories across individuals (12). Disregulated immune activation or suppression may lead to chronic inflammation, disruption of immune rhythms, and even immunopathological damage. The bidirectional regulation between the pro-inflammatory cytokine IL-6 and immunosuppressive pathways is central to limiting excessive inflammation and maintaining homeostasis, yet under certain conditions it can also drive the system toward a pathological state of sustained high inflammation (13, 14). Thus, it is necessary to modeling the immune network as an integrated system, particularly with humoral, cellular immunity and IL-6 interacts (15). Existing studies of within-host viral dynamics are generally based on the target cell limited framework, which has established a standard modeling paradigm for describing viral infection, replication, and clearance of susceptible cells, and has enabled the estimation of key parameters through data fitting (16). In the context of acute infections, innate immunity is often incorporated to explain the decline and clearance phases of viral load (17). Additional extensions encompass adaptive immunity, wherein cellular (CTLs/NK cells) and humoral (antibodies) components collectively facilitate targeted clearance and neutralization (15). More comprehensive models may also include quantitative analyses of concurrent innate and adaptive immunity actions linking infection dynamics to pathological outcomes (18, 19). However, most studies used models with only one or two immune pathways (20, 21), leaving underexplored the interactions between pro-inflammatory cytokines and suppressive mechanisms, as well as scenarios involving exogenous viral inputs.

In this study, a virus-immune interaction network model is proposed to integrate viral dynamics with key host immune processes including virus [V], innate immunity [I], cellular immunity [C], humoral immunity [H], immune suppression [S], and pro-inflammatory IL-6 cytokine [IL6]. As this study highlights the qualitative dynamical behavior, all state variables represent the dimensionless strength rather than the concrete quantitative numbers. So the viral load, immune cell populations, and IL-6 levels thus represent scaled relative values rather than absolute concentrations. All quantities in equations and figures appear in dimensionless units. Hill functions are employed for unified description of cross-module activation, inhibition, and saturation effects (22). This approach captures concentration-dependent viral replication and nonlinear inhibition by multiple immune clearance pathways. In addition, an explicit exogenous viral input function α(t) is introduced to simulate different infection scenarios, including persistent and transient exposures. This design preserves the biological interpretability of both variables and parameters. Based on the proposed virus–immune interaction network model, we systematically investigate the effects of external viral input patterns, key process parameters, and immune module combinations on the overall dynamical behavior. We first examine system responses under conditions of persistent viral input, showing that different combinations of immune clearance efficiency and viral replication capacity may lead to steady-state transitions, multistability, or the emergence of new pathological states. Critical parameter thresholds and their associated steady-state changes are quantified through numerical simulations and bifurcation analysis. We then analyze the recovery dynamics under non-persistent viral input scenarios, and find that different immune modules exhibit distinct timescales during the post-clearance recovery process, with lagged declines observed across modules. In addition, to quantitatively describe the temporal features of the infection and inflammation phases, we introduce two time-based indicators, “infectious duration” and “illness duration” and assess their sensitivity to immune regulatory mechanisms under multiparameter perturbations. Finally, through subsystem enumeration analysis, we systematically screened variable combinations while fixing partial modules at steady state, identifying the triplet of virus ([V]), innate immunity ([I]), and cellular immunity ([C]) as the core structure driving complex dynamical modes such as bistability and oscillations. In contrast, humoral immunity and the suppression–inflammation loop primarily regulate response amplitude and recovery synchrony.

In summary, this work reveals the mechanisms underlying multistability and oscillations in the immune network, while providing quantifiable temporal indicators for evaluating the efficiency of immune regulation. It advances the theoretical understanding of virus–immune system coupling, offerring a scalable analytical framework and methodological basis for the design and parameter optimization of personalized immune intervention strategies.

2 Methods

2.1 Model discription

Significant differences exist in immune responses among individuals following viral infection, and these differences are closely associated with clinical manifestations. To elucidate the pathogenic mechanisms of viral infection (23, 24) and to provide theoretical insights for therapeutic strategies (25), we developed a virus–immune interaction model based on ordinary differential equations, with a particular focus on the overall regulatory role of the immune system. The model adopts a modular structure to represent the interactions among virus, innate immunity, adaptive immunity, immune suppression, and inflammatory cytokines, aiming to highlight the nonlinear couplings across multiple feedback loops. The present model is formulated as a coarse-grained representation of a highly complex immune regulatory network. Accordingly, each dynamical variable reflects the integrated intensity of the functional response of a corresponding immune module, rather than the exact number or concentration of a specific cell subtype or molecular species. Steady-state and bifurcation analyses are employed to reveal the possible dynamical behaviors of the system. In contrast to recent approaches using spatially explicit or hybrid dynamical models to explore immune spatial features (26), this study emphasizes module coupling and stability analysis at the global level.

The human immune system contains roughly two components, innate immunity and adaptive immunity (1, 27). Upon viral invasion, innate immunity serves as the first line of defense, rapidly eliminating pathogens through physical and chemical barriers (e.g., skin, gastric acid), phagocytic cells (e.g., macrophages and neutrophils), and natural killer (NK) cells. This response is accompanied by an increase in pro-inflammatory factors such as IL-6, which establishes the inflammatory and signaling background for subsequent responses. Adaptive immunity is then activated to achieve targeted clearance, providing long-lasting protection against specific antigens (28). Within adaptive immunity, T cells play a central role that CD8⁺ T cells recognize and destroy host cells infected by viruses or other pathogens, while CD4⁺ T cells coordinate interactions among immune cells and assist in their activation. Humoral immunity relies on B cells and their differentiation into plasma cells, which produce antibodies to neutralize pathogens and prevent further damage. At the same time, to avoid excessive immune activation, suppressive pathways provide negative feedback on inflammation. Regulatory T cells (Tregs) inhibit overactive T cells, while immune checkpoint proteins (e.g., PD-1, CTLA-4) and anti-inflammatory cytokines (e.g., IL-10, TGF-β) also play critical roles in dampening interactions among immune cells. In the present model, these mechanisms are incorporated as interactions among the corresponding immune modules.

In the regulation of the immune system, humoral immunity also plays an important role. Humoral immunity influences the host’s defense capacity against pathogens by modulating immune cell functions and the strength of interactions among immune cells. Hormones such as cortisol, sex hormones (estrogen and testosterone), and thyroid hormones can regulate immune cell interactions by either promoting or suppressing immune cell activity and differentiation. Cortisol exerts a significant immunosuppressive effect by attenuating the strength of immune cell interactions, whereas estrogen generally enhances such interactions, particularly by promoting antibody production through B cells (29, 30). The bidirectional regulation of IL-6 and suppressive mechanisms plays a key role in limiting excessive inflammation, maintaining homeostasis, and coordinating temporal responses. IL-6 is a major pro-inflammatory cytokine (31) that promotes immune cell activation and differentiation, especially in antibody production and T-cell function. It also plays a central role in the acute-phase response by stimulating hepatic synthesis of C-reactive protein (CRP) and enhancing local inflammatory reactions. Elevated IL-6 levels are closely associated with the development and progression of various diseases, including infections, cancers, and autoimmune disorders. Notably, IL-6 not only drives inflammatory responses but can also induce the activation of immunosuppressive pathways, thereby establishing a dynamic balance between pro-inflammatory and suppressive effects. This dual role makes it a critical node in regulating infection outcomes and immune homeostasis.

To effectively characterize the host immune response mechanisms following viral invasion, we construct an interaction network of viral infection and immune response comprising six variables that represent different immune modules (Figure 1). This network accounts for the complex interactions among viral dynamics ([V]), innate immunity ([I]), cellular immunity ([C]), humoral immunity ([H]), immune suppression ([S]), and interleukin-6 ([IL6]) after infection, thereby capturing the global dynamical features of the host immune response following viral entry.

Figure 1

In the model, the dynamic evolution of viral load results from multiple interacting mechanisms. First, the virus may enter the host system through external sources. To capture this effect, a time-dependent input function, α(t), is introduced. During a finite initial period [0, T], the function is maintained at a constant value ϕ, representing continuous exogenous viral exposure; thereafter, it drops to zero, indicating that external infection sources cease to contribute. In this way, α(t) characterizes the environmental exposure faced by the host at different stages. The mathematical definition of α(t) is provided in Equation 1.

In addition to exogenous input, the virus itself is capable of replication. The proliferation rate of the virus does not increase linearly but instead exhibits a nonlinear saturating behavior as the concentration rises. To capture this concentration dependence, the model employs a Hill function (Equation 2) to characterize the kinetics of viral growth.

Here Equation 2 reflects the trend that the viral replication rate increases with viral abundance, revealing the saturating nature constrained by resources, environmental factors, or host conditions (22). a_v0 denotes the maximum replication rate of the virus, i.e., the limiting speed under saturation; n_v0 is the Hill coefficient, which determines the steepness of the process and reflects the sensitivity of replication rate to changes in viral concentration; and k_v0 is the half-saturation constant, corresponding to the viral concentration at which the replication rate reaches half of its maximum value. In other words, this term represents the kinetic approximation of the viral self-replication processes. Simultaneously, the host immune system continually participates in viral clearance, as innate immunity, cellular immunity, and humoral immunity all accelerate viral elimination. Specifically, virus death induced by cellular immunity is described by Equation 3.

Similarly, viral clearance mediated by innate immunity and humoral immunity can be represented by Equation 4.

In Equations 3 and 4, each term corresponds to the clearance effect of a specific immune mechanism, dv₁, dv₂, and dv₃ respectively quantify the clearance efficiencies of cellular immunity, innate immunity, and humoral immunity against the virus. The associated n and k parameters characterize the sensitivity and half-saturation level of these immune effects. The multiplication by [V] represents the dependence of the viral clearance rate depends on the strength of immune responses and on the viral load itself. Finally, even in the absence of immune action or exogenous input, viruses gradually decline due to intrinsic inactivation, environmental instability, and passive removal. This natural decay effect is modeled by d_v4[V], serving as a background dissipative mechanism. The complete expression for viral dynamics, as shown in Equation 5, is obtained by integrating these four components of exogenous input, self-replication, immune clearance, and natural decay.

Innate immunity ([I]) describes the rapid defensive response of the host during the early stage of viral invasion. Its dynamical evolution integrates four processes, including viral activation, self-amplification of immunity, negative feedback regulation by suppressive mechanisms, and natural decay of immune factors. First, the presence of the virus directly triggers the activation of innate immunity, which exhibits a saturation effect. This effect is represented using a concentration-dependent Hill function, as shown in Equation 6.

where, denotes the maximum activation rate triggered by viral stimulation; is the Hill coefficient, characterizing the steepness of the innate immune response to changes in viral level; and is the half-saturation constant, corresponding to the viral load required for immune activation to reach half of its maximum level. This process reflects the sensitivity of host immune cells or molecules to viral stimulation during the early stage of infection. Second, the innate immune system possesses a certain degree of self-amplification. For example, interferon molecules, once activated, can further enhance their own production, thereby forming a positive feedback loop. This effect can be represented by Equation 7.

where represents the maximum self-activation rate, while and respectively determine the sensitivity of the response and the half-saturation point. This term reflects that, once initiated, innate immunity can be rapidly amplified through positive feedback, thereby strengthening its ability to suppress the virus. At the same time, the immune system must also avoid excessive reactions. To this end, the model introduces an immune suppression module ([S]), which suppresses innate immunity via negative feedback. This regulation can be written as follow.

where represents the strength of the suppressive effect, while and control the nonlinear threshold characteristics of the process. Through this mechanism, the model reflects that host defense, while being activated, is also subject to regulation in order to prevent excessive immune responses that could cause tissue damage. Finally, innate immune factors undergo natural death or degradation, represented by the term , which indicates that immune molecules or cells gradually decay in the absence of continuous stimulation. By combining Equations 6–8, the dynamics of innate immunity can be summarized in the form of Equation 9. This equation comprehensively characterizes the activation, amplification, regulation, and decay mechanisms of innate immunity, reflecting the rapid and dynamic defensive features exhibited by the host during the early stage of infection.

During viral infection, cellular immunity ([C]) mainly refers to the specialized defense provided by T cells. Its core function is to recognize and eliminate host cells infected by the virus, thereby blocking viral replication and transmission (32). The dynamics of this process can be decomposed into four aspects, namely initiation, amplification, regulation, and death. First, the activation of cellular immunity requires a “dual condition”, namely the presence of a sufficient viral load ([V]) together with the prior activation of innate immunity ([I]). Accordingly, Equation 10 effectively characterizes the early activation process of cellular immunity.

where denotes the maximum initiation rate, while and respectively control the sensitivity and threshold of cellular immunity to the levels of innate immunity and viral concentration. The multiplicative form of this term clearly reflects the fine regulation of the immune system, ensuring that a specific cellular immune response is effectively initiated only when the viral concentration exceeds a certain level and innate immunity has been activated. Second, once cellular immunity is initiated, a positive feedback mechanism is also present (33). Activated T cells or effector cells can enhance their own expansion and maintenance through cytokine-mediated processes. This effect is represented in the form of a Hill function, as shown in Equation 11.

where represents the maximum self-activation rate, and respectively determine the sensitivity and half-saturation level of the response. This term reflects that, once initiated, cellular immunity can form a sustained amplification effect, maintaining a high level of viral clearance ability over an extended period. Meanwhile, cellular immunity is also regulated by immune suppressive factors ([S]). The negative feedback mechanism governing this regulation is represented by Equation 12.

where represents the suppression intensity, and and characterize the nonlinear features of this process. This term reflects the limitation of cellular immunity activity by immune suppressive factors, preventing excessive immune responses that could cause immune-related tissue damage. Finally, effector cells of cellular immunity also undergo natural apoptosis or inactivation, a process represented by the term , indicating that even in the absence of suppressive factors, effector cells such as T cells gradually decline. Therefore, by integrating the above Equations 10, 11, and 12, along with the natural apoptosis of cellular immunity, the dynamics of cellular immunity in viral infection can be represented by Equation 13.

This equation systematically encapsulates the entire process of cellular immunity during viral infection, driven by both innate immunity and the virus, including activation, self-amplification, negative feedback regulation by suppressive mechanisms, and eventual natural decay.

In the host immune response, humoral immunity ([H]) primarily represents the adaptive immune response mediated by B cells. Compared to innate immunity and cellular immunity, the activation of humoral immunity typically occurs later during infection; however, it can effectively clear free viral particles through the specific recognition and neutralization by antibodies (34). Its dynamics integrate activation signals from cellular and innate immunity, are regulated by immune suppression mechanisms, and are accompanied by the natural decay of humoral immunity itself. One major source of activation for humoral immunity is cellular immunity ([C]). T cells directly kill infected cells and promote B cell differentiation and antibody production through cytokine secretion. This process is represented by Equation 14.

where represents the maximum facilitation rate of cellular immunity in activating humoral immunity, while and respectively determine the nonlinear steepness and half-saturation point of the response. This term reflects the synergistic relationship between cellular and humoral immunity, where T cells provide critical support for B cell activation and antibody production. Furthermore, innate immunity ([I]) can also provide activation signals for humoral immunity. For example, interferon molecules or inflammatory factors can indirectly enhance B cell function. This effect is represented by Equation 15.

here, represents the maximum facilitation rate of innate immunity in activating humoral immunity, while and control the dynamics of this process. This term reflects the cascade relationship between different immune modules, enabling humoral immunity to be effectively activated under the influence of multiple signals. Meanwhile, humoral immunity is also regulated by negative feedback from the immune suppression module ([S]), as shown in Equation 16.

where represents the suppression intensity parameter, while and determine the threshold characteristics of this effect. This term reflects that, while activating humoral immunity, the immune system also prevents excessive antibody responses through suppressive mechanisms, thereby avoiding unnecessary damage to the host. Finally, humoral immune molecules (such as antibodies) undergo a natural degradation process, represented by the term , which indicates that, even in the absence of viral stimulation, antibody levels gradually decline. Therefore, the dynamics of humoral immunity can be summarized in the form of Equation 17, which comprehensively reflects the key role of humoral immunity in adaptive immune responses. Specifically, it is initiated by the combined activation signals from cellular and innate immunity, exerts neutralizing effects of specific antibodies under the regulation of immune suppression, and gradually decays after the response ends.

In the immune system, it is not possible for all modules to maintain a high state continuously. Therefore, immune suppression regulates this by introducing a negative feedback mechanism through the immune suppression module ([S]), which reflects the host’s effort to prevent overactivation of immune responses. This module plays a crucial role in maintaining immune homeostasis and preventing tissue damage. The dynamics of this process are influenced by activation and regulation from cellular immunity ([C]), innate immunity ([I]), and inflammatory factor IL-6, while also being limited by its own natural decay. First, immune suppression has a baseline level, represented by the constant δ. In the model, δ = 0.0001, which signifies that even in the absence of external stimuli, the host maintains a certain degree of immune suppression to ensure homeostasis. Second, under external stimulation, both cellular immunity ([C]) and innate immunity ([I]) can promote the production of immune suppressive factors. This process is represented by Equation 18.

where and represent the maximum rates at which cellular immunity and innate immunity activate immune suppressive factors, respectively, while , and , control the nonlinear characteristics of their responses. This part reflects that when the immune system is active, the body simultaneously enhances immune suppression to avoid excessive immune responses. In addition, the inflammatory factor IL-6, at high levels, further drives the suppressive effect, forming an inflammation suppression negative feedback loop, as shown in Equation 19.

where represents the intensity parameter for this feedback effect, while and characterize its threshold and nonlinear effects. This term reflects that when inflammation increases, the body attempts to enhance the suppressive response in order to reduce inflammatory damage. Finally, immune suppressive factors undergo natural decay or inactivation, represented by . This term indicates that, even in the absence of continuous activation, the level of immune suppression gradually declines.

In summary, the dynamics of the immune suppression module can be expressed by Equation 20. This equation incorporates both the baseline level and the upregulation effect from immune activation, while also considering the negative feedback regulation driven by inflammatory factors and the natural decay of suppressive factors, thereby maintaining a dynamic balance between activation and suppression in the immune system.

Interleukin-6 (IL-6) is an important pro-inflammatory cytokine produced during viral infection, widely regarded as a key marker of immune system activation and inflammation levels (31). The dynamic changes in its levels reflect the intensity of the host’s immune response to infection and inflammation-associated pathological conditions. In the model, the dynamics of IL-6 are primarily regulated by viral load, cellular immunity levels, and its own degradation process. First, cellular immunity ([C]) is one of the key drivers of IL-6 production. Activated T cells and associated immune factors promote the secretion of IL-6, a process represented by Equation 21.

where represents the maximum activation rate, while and describe the sensitivity and half-saturation level of this process, respectively. This term reflects the pro-inflammatory effect associated with active cellular immunity. Additionally, viral load ([V]) itself can directly stimulate the production of IL-6, especially during infection spread and inflammation escalation. This process is described by Equation 22.

where represents the maximum rate of IL-6 production induced by the virus, while and determine the nonlinear characteristics of the response. This term reflects the positive correlation between viral load and inflammation levels. Finally, IL-6 molecules undergo a natural degradation process, represented by . Even in the absence of stimulation, IL-6 gradually declines to maintain homeostasis. Therefore, the dynamics of IL-6 can be uniformly expressed by Equation 23.

This equation captures the role of IL-6 in infection and inflammation, where its production is driven by cellular immunity and viral load, and its accumulation is limited by self-degradation. A sustained increase in IL-6 levels typically indicates the onset of a high-inflammation state in the host, potentially leading to severe pathology and clinical symptoms.

The model equations are summarized in Equations 24.

2.2 Parameters and numerical solution

Firstly, the well-posedness of the model (24). The right-hand side of the system consists of linear terms, bilinear interactions, and Hill-type nonlinearities of the form , all of which are continuously differentiable and therefore locally Lipschitz on . By the classical existence and uniqueness theorem for ordinary differential equations, these properties ensure that the system admits the existence of a unique solution for the biologically meaningful initial condition problem. Furthermore, the boundedness and positivity of the solution is ensured by the vector field of the ODEs because the decay or clearance terms are proportional to the corresponding variable and vanish at zero such that the vector fields always points to the positive direction, preventing trajectories from leaving the biologically relevant domain.

To better capture the heterogeneity of host immune responses during viral infection, both in terms of intensity and timing, and to account for the resulting diversity of clinical manifestations, we systematically explored the parameter space of the virus–immune response interaction network and performed corresponding dynamical simulations. To balance exploration capacity with computational efficiency, parameter-space sampling was implemented. The model equations include six variables and 63 parameters, making it impractical to fully cover the high-dimensional parameter space; therefore, the sampling dimensions and sample size were reduced to improve computational efficiency. Specifically, the half-saturation constants and Hill coefficients, which characterize system dynamical features, were fixed and not sampled. In bifurcation analyses, random sampling was applied to major parameters, including the viral replication rate, activation and clearance rates among immune compartments, decay rates of immune compartments, and cytokine production rates. For each parameter set, initial values of the variables were assigned, and the system of ordinary differential equations was numerically solved using the ode15s solver in MATLAB.

Furthermore, we examined how different viral input patterns influence the dynamical trajectories of immune state variables. To investigate how qualitative behaviors change under parameter variation, we also performed bifurcation analysis using MatCont, which enabled us to detect stability shifts and regions where multiple steady states coexist. In addition, we assessed the dependence of long-term outcomes on initial immune-state values and used kernel density estimation to visualize the distribution of steady states reached under different initial conditions.

In addition, to ensure that the model computations capture the qualitative characteristics of the system, appropriate initial conditions and parameter values were selected based on existing immunological studies. Specifically, the initial conditions specify the concentrations of virus, immune cells, and cytokines set at the beginning of the simulations. These concentrations represent the baseline state of the immune system, i.e., the normal levels in the absence of viral infection or other external perturbations.

For this purpose, all state variables were treated in a dimensionless form prior to simulation. The present model is constructed as a coarse-grained representation of a highly complex immune regulatory network. Each state variable does not correspond to the exact number or concentration of a specific cell type or molecular species, but rather represents the integrated intensity of the functional response of the corresponding immune module. In real physiological processes, each immune function typically involves the cooperation of multiple cell subtypes and signaling pathways in different activation states, whose overall effect cannot be uniquely and unambiguously mapped to a single measurable cell numbers or molecular levels. Given that this study focuses on qualitative dynamical analysis and aims to explore the interaction structure among immune modules and the resulting system-level behaviors, rather than to perform quantitative data fitting, the modeling objective is to capture the essential qualitative features of the system. For the present framework, retaining dimensionless state variables is consistent with its conceptual objective and avoids implying a level of quantitative precision that the model is not designed to provide. Accordingly, only the time scale is specified, the units for other state variables are dimensionless while no explicit physical units are assigned to the state variables. Based on this dimensionless quantity, the parameter values involved in this study should be interpreted as rescaled quantities in the normalized model, rather than as directly measurable physiological rates or concentrations. Specifically, while the time variable is kept in physical units (days), as stated above, the state variables are dimensionless, and the remaining coefficients represent aggregated interaction strengths and effective clearance/activation rates at the module level. Therefore, this parameter set is intended to reproduce qualitative dynamical mechanisms, rather than to provide quantitative fitting of unit-consistent estimates.

In contrast, the parameter values are used to characterize the key mechanisms of immune responses, including how the virus infects host cells and how immune cells are activated or suppressed. These values were chosen according to the known characteristics of immune responses. For example, these include the sequential order of immune module activation following viral invasion, the strength of interactions among immune cells, the modes of action of suppressive cytokines, and the mechanisms by which cytokines regulate immune responses (1, 35–37). The specific parameter settings are listed in Table 1.

3 Results

3.1 Characteristic dynamical behavior

Based on the above setup, to elucidate the dynamical interplay between viral replication and immune responses, we performed numerical simulations to examine system evolution under different conditions and selected representative outcomes of virus–immune interactions. Figures 2a1, a2 illustrates the dynamics of immune modules under continuous external viral input. During the early phase of sustained exposure, immune modules are rapidly activated, while the suppression module decreases, collectively resisting viral invasion. Although the immune modules can transiently restrain viral growth at first, viral expansion arises from external input as well as self-replication, which is inherently nonlinear. The half-saturation property of the Hill function limits unbounded viral proliferation, yet under persistent external input, even a small influx may push the system beyond a latent threshold. Once this threshold is crossed, the existing immune modules become insufficient to control viral growth, leading to a rapid rise in viral load. In response, immune modules adjust accordingly, and the system eventually converges to a new steady state. This phenomenon suggests that the host exhibits a latent period under sustained viral exposure, after which viral load escalates rapidly. These results indicate that external viral input can trigger an irreversible transition of immune system states, reflecting the system’s multistability.

Figure 2

Under short-term external input, the system returns to its original steady state after the input is removed (see Figure 2b1), indicating that a brief viral exposure can induce transient fluctuations in the immune system, but the overall structure remains stable without state transition. In contrast, when the input duration is sufficiently long (e.g., applied for 70, units before removal), the system trajectory is pushed into another basin of attraction and eventually converges to a new steady state. Subsequently, when viral input is applied again during t ∈ [140, 210], the input no longer alters the system’s steady state (see Figure 2). This demonstrates that although viral input is removed after a period of time, the accumulated effect has already driven the system to a new steady state, such that later removal or reintroduction of input produces no significant change in the system’s stability.

A more special case is shown in Figure 2c, which illustrates that under certain initial conditions and parameter configurations, the system may enter a sustained oscillatory state. In this case, variables such as virus load [V], innate immunity [I], cellular immunity [C], and humoral immunity [H] all exhibit stable periodic oscillations. This dynamical pattern typically indicates that the immune system, after viral input, is driven into a critical regime, where the interplay of positive and negative feedback among different immune factors leads the system into periodic fluctuations. It reflects both the sensitivity of the system to viral input and its inherent nonlinear dynamics, and may correspond to pathological phenomena such as chronic inflammation or immune rhythm disorder.

3.2 Bifurcation behavior of the system

To characterize how the immune–virus system transitions between different qualitative regimes, we first outline the stability criterion used throughout the bifurcation analysis. Let the viral–immune dynamics be written in vector form as = F(x), x = (V, I, C, H, S, L)^T. An equilibrium point satisfies . To analyze its local behavior, we linearize the system around as , where the Jacobian matrix is defined by . According to the linearization principle and the Hartman-Grobman theorem, the eigenvalues {λ_i} of J(x) fully determine the stability of the equilibrium, if Re(λ_i)< 0 for all i, the equilibrium is locally asymptotically stable; if at least one eigenvalue satisfies Re(λ_i) > 0, the equilibrium becomes unstable; if an eigenvalue crosses zero, the system undergoes a saddle-node or transcritical bifurcation; and if a complex conjugate pair crosses the imaginary axis, i.e., λ₁,₂ = ± iω, then a Hopf bifurcation occurs, leading to the emergence of periodic oscillations (38). These stability properties are particularly relevant for immune–virus dynamics, where nonlinear activation and suppression pathways can generate several coexisting steady states. Transitions between these states occur when the leading eigenvalue of the Jacobian changes sign, while oscillatory responses emerge as feedback loops bring a conjugate pair of eigenvalues close to the imaginary axis, eventually giving rise to Hopf bifurcations. In this work, the Jacobian is evaluated along each solution branch to assess stability and to identify points where such transitions occur. Branch segments with all eigenvalues having negative real parts are plotted as solid curves, whereas segments containing at least one eigenvalue with positive real part are shown as dashed curves (39). This provides a coherent framework for interpreting the multistability, switching behaviors, and oscillatory regimes that appear in the subsequent bifurcation diagrams. With this framework, we next examine how the system responds under continuous viral input and how parameter variations shape the qualitative structure of the immune dynamics.

To analyze the dynamical characteristics of the system under continuous external viral input, we first consider the case of α(t) = ϕ. The results show that, with parameters fixed, the system can converge to different steady states solely by varying the initial conditions. As shown in Figures 3a1, a2, panel (a1) depicts the time trajectories of viral load V(t) under multiple initial conditions, while panel (a2) shows the trajectories of cellular immunity C(t). The black and red dashed lines indicate the times when the system reaches steady states, t₁ = 7 and t₂ = 13, respectively, demonstrating a negative correlation between the response time and the response intensity of cellular immunity (40). It can also be seen that although variable innate immunity [I] ultimately converges to the same steady state, cellular immunity [C] reaches two distinct stable states. Furthermore, under continuous viral input (ϕ = 0.01), to examine the influence of key parameters on the system’s final state, we performed a bifurcation analysis of the model while keeping the initial conditions fixed. Figures 3b1 – b6 display the branch structures of stable and unstable solutions as key parameters vary. Different colors represent different steady-state solution branches, with solid lines denoting stable solutions and dashed lines denoting unstable solutions, clearly revealing the existence of multistability regions and their dynamic pathways.

Figure 3

It is particularly noteworthy that the parameter d_v1 in the model equation specifically represents the efficiency of cellular immunity in clearing the virus, reflecting the key impact of changes in the effectiveness of the cellular immune mechanism on the viral infection system. From the results shown in the figures, it is clearly observed that as the cellular immunity clearance ability (d_v1) changes, the system exhibits significant multistability and bifurcation behavior. Figures 3b1–b6 illustrates the multistability and bifurcation features of the system regulated by d_v1, underscoring the pivotal role of cellular immunity in shaping global immune steady states. Specifically, when the clearance ability of cellular immunity is weak (small d_v1), the viral steady-state level remains high, accompanied by elevated levels of the inflammatory cytokine IL-6, corresponding to clinical scenarios such as chronic viral infection or impaired host immune function. As d_v1 gradually increases, the enhanced clearance efficiency of cellular immunity drives a marked reduction in viral steady-state levels. Crossing one or more critical bifurcation points, the system transitions from a high-viral steady state to a low-viral steady state, entering a non-infectious healthy state. This process reflects a recovery pattern in which enhanced immune capacity or therapeutic intervention progressively clears the virus. Furthermore, near the bifurcation points, the system is highly sensitive to initial conditions, that is, different initial immune levels can lead to drastically different disease progression trajectories, manifested as either a long-term chronic infection state or a rapid recovery state with viral clearance. Specifically, under the same parameter conditions, both high viral-high inflammation steady states and low viral-low inflammation steady states coexist. At this point, the host’s initial immune level (e.g., initial immune cell count, baseline inflammatory response) will determine which steady state the system ultimately converges to. This phenomenon suggests that even with the same treatment, individuals may experience entirely different disease outcomes due to variations in their initial immune conditions.

To further elucidate the joint regulatory effects of viral replication and clearance mechanisms on the global stability structure of the system, this study extends the single-parameter bifurcation analysis by selecting two key dynamical factors including the viral replication rate parameter and the viral clearance rate parameter . A two dimensional parameter space scan was conducted to systematically evaluate their influence on the multistability structure of the model (41) (see Figure 3c).

Specifically, all other parameters were fixed, and a wide-range scan of and was performed in double logarithmic coordinates. For each parameter combination, the number of steady states and their stability properties were determined by solving the steady-state system and performing eigenvalue analysis. The colors in the figure represent the number of steady states, ranging from monostability (blue) to five steady states (orange-red), showing a distinct regional distribution as illustrated in Figure 3c. It can be observed that in the region where is small and is large, the system possesses only one stable solution, corresponding to a healthy state in which the virus is effectively controlled by the immune system. When the viral replication rate increases (larger a_v0) and the clearance efficiency decreases (smaller d_v1), multiple steady states emerge, indicating that the system dynamics may fall into different infection outcomes, including viral persistence and uncontrolled immune activation. In the transitional regions, rich tri-stability and multistability structures are observed, suggesting strong sensitivity to initial conditions and nonlinear response characteristics. In particular, in the orange regions, the system exhibits five steady states, reflecting highly complex steady-state topologies that may be accompanied by abundant bifurcation phenomena and unpredictable dynamical responses. This two-parameter scan reveals the global steady-state distribution of the system under different combinations of viral replication and clearance intensities, emphasizing the nonlinear complexity of virus–host immune interactions. The results suggest that, within certain parameter ranges, the system may exhibit multiple potential courses of infection, where even small differences in the host’s initial immune state or intervention strategies can lead to completely different clinical outcomes. Therefore, this finding highlights the importance of individualized prediction and treatment of infectious diseases.

3.3 Coexistence of higher and lower steady states

In Figure 3c, the steady states change very rapidly for certain parameter values. When d_v1 is large (e.g., d_v1 = 100), multiple changes in steady states occur as k_v0 increases. In fact, immune responses during viral infection vary significantly among individuals and are influenced by factors such as age, health status, and sex (42, 43). To further investigate the stochasticity and robustness of the system, several representative parameter points in specific regions were selected for analysis. As shown in Figures 4a–c, when α(t) is continuously applied and d_v1 is large, while keeping all other parameters unchanged, the temporal dynamics of viral load [V], innate immunity [I], cellular immunity [C], and interleukin-6 (IL-6) were examined. The results indicate that, due to the relatively large value of the viral natural degradation coefficient d_v4, the viral load rapidly decreases to very low levels, or even close to zero, across all initial conditions. This outcome is consistent with the biological expectation of rapid viral clearance. However, other immune modules ([I], [C], [IL6]) exhibit different trajectories. In particular, within the cellular immunity module, one trajectory remains at a high state, while another stays near zero, corresponding to a low state. We refer to these two scenarios as the high state and low state of cellular immunity (44).

Figure 4

When cellular immunity is in the high state, immune responses (such as T-cell activity) serve as a direct mechanism against viral infection. In this case, cellular immunity is highly active, allowing the immune system to effectively recognize and eliminate virus-infected host cells through T cells or other effector cells. Although the viral load is V = 0, cellular immunity [C] and the inflammatory cytokine IL-6 remain in the high state. This may reflect a “lag effect” of the immune system (45) or the self-sustaining nature of immune responses. During viral clearance, the high states of cellular immunity [C] and IL-6 may not be entirely virus-driven but instead maintained through their own positive feedback loops. Sustained high expression of cellular immunity [C] can further activate and maintain innate immunity [I] at relatively high levels, forming a positive feedback cycle of immune activation. In addition, as a key inflammatory mediator, IL-6 is directly or indirectly activated by cellular immunity [C], and through positive feedback further enhances immune cell activity. Such positive feedback mechanisms may trap the system in a prolonged high-expression state. Thus, even when the viral load is zero, an excessively activated immune system can still induce inflammation. These high-expression states typically occur shortly after infection and represent part of the acute-phase immune response. Without external intervention, patients may remain in chronic pathological conditions, such as persistent inflammation or immune overactivation. In contrast, when cellular immunity is in the low state, the virus is still effectively cleared, but cellular immunity [C] and the inflammatory cytokine IL-6 rapidly decline to near-zero levels. This may reflect the action of immune suppression mechanisms or other immune modules (such as humoral and innate immunity) that effectively shut down immune responses after viral clearance, thereby preventing tissue damage caused by excessive immune activity. At the same time, innate immunity [I] and IL-6 also decrease rapidly when cellular immunity [C] is low, showing coordinated reduction. This indicates that strong cooperative regulatory mechanisms exist in the system, where the low activity of cellular immunity is further reinforced by negative feedback or immune suppression pathways, stabilizing the shutdown of immune responses. This state represents a more favorable recovery of immune homeostasis, suggesting that the host can return to normal levels quickly and reduce the risk of chronic inflammation or immune-related disorders. In summary, the system exhibits a clear bistable behavior. This phenomenon reflects the diversity and complexity of interactions among immune cells and further demonstrates that the ultimate steady state of the immune system is strongly influenced by the initial immune conditions. These findings highlight the importance of early immune regulation in preventing long-term inflammation or chronic immune activation and provide theoretical support for the development of early intervention strategies in clinical practice.

To further analyze the sensitivity of immune system steady states to initial conditions, a state classification diagram was constructed with the initial level of innate immunity I₀ on the horizontal axis and the initial level of cellular immunity on the vertical axis, systematically examining the influence of different initial configurations on steady-state outcomes (Figure 4). Here, the viral clearance efficiency of cellular immunity was fixed at a relatively high level to ensure that the virus could be cleared within a short time, thereby eliminating the interference of continuous viral input on immune steady states.

In the figure, red, blue, and green represent three typical combinations of steady-state responses. Red corresponds to the case where cellular immunity [C] and the inflammatory cytokine IL-6 are both in high-expression states, indicating the possibility of immune overactivation or chronic inflammation. Blue denotes elevated IL-6 while cellular immunity [C] remains in a low state, suggesting that inflammation is activated but cellular immunity is not effectively engaged. Green represents both variables at low-expression levels, corresponding to the ideal state in which the virus is completely cleared and the immune response is successfully shut down.

Figures 4d–e show the discrete distribution of final steady states under different combinations of initial values. The three types of steady states exhibit clearly defined boundaries in the initial state space (black dashed lines). However, in certain boundary regions, noticeable overlaps are observed, indicating that different steady states may arise from similar initial conditions. This phenomenon reflects the multistability characteristics of the system.

To further characterize the distribution of different steady states in the initial condition space, we computed the occurrence frequency of each steady state on the parameter plane and presented the results as a density distribution map, as shown in Figure 5. As a complement to the discrete classification results in Figures 4d, e, this figure provides a more intuitive reflection of the probabilistic distribution of the three steady states in the boundary regions. The results show that while different steady states are clearly separated in most regions, significant overlaps and intersections occur near the critical boundaries. This indicates that the final outcome of the system is highly uncertain under such initial conditions, exhibiting typical features of multistability and sensitivity to initial conditions. In other words, the density distribution map provides a continuous probabilistic perspective for the discrete results, further highlighting the complex dynamical features in which the system may evolve into different steady states even under identical initial conditions. The results indicate that, when exposed to the same viral input, the long-term steady state of the immune system is sensitive to both system parameters and initial conditions. Particularly, during the early stages of disease progression, only small differences in the immune system state may lead to markedly different inflammatory outcomes.

Figure 5

3.4 Immune recovery dynamics and analysis of duration indicators under finite viral input

To investigate the recovery dynamics of the immune system under finite viral input, we simulated the case in which the virus was continuously introduced for 7 days and then stopped. as is shown by Equation 25. As shown in Figure 1b1, ϕ = 0.01, the viral load continuously increases and reaches a peak before the termination point of input at t = 7. After the input ceases, the viral load rapidly declines and approaches zero, indicating that the virus is gradually cleared by the immune system.

In terms of the responses among different immune modules, innate immunity ([I]), cellular immunity ([C]), and humoral immunity ([H]) are all activated during the viral input period, showing a clear upward trend. Among them, humoral immunity ([H]) remains at a high level even after viral input ceases, then gradually declines and stabilizes, exhibiting a typical delayed-response characteristic. As a representative marker of the inflammatory response, IL-6 rises synchronously during the viral peak but decreases much faster than the immune cell indicators. This phenomenon suggests that the inflammatory response decays rapidly during viral clearance, whereas the decline of immune cell activity shows a pronounced lag (46).

From the perspective of the overall system dynamics, the immune variables gradually return to steady states about 5 days after the termination of viral input (t = 12), recovering to levels close to the initial immune state. This indicates that the immune system possesses strong self-regulatory capacity. However, the steady-state values of humoral immunity ([H]) and the immune suppression module ([S]) remain slightly higher than their initial levels and do not fully return.

This phenomenon can be explained from both the model mechanism and immunophysiological perspectives. From the viewpoint of dynamical structure, humoral immunity [H] is jointly activated by cellular immunity ([C]) and innate immunity ([I]). After viral clearance, the residual activities of cellular immunity [C] and innate immunity [I] can still drive humoral immunity humoral immunity [H] to remain at a relatively high level for a short period. Together with the relatively small decay rate of humoral immunity [H], this results in a delayed decline. Similarly, the production of immune suppression [S] is driven by a constant term and immune activation variables (cellular immunity [C] and innate immunity [I]), while being suppressed by the inflammatory factor IL-6. After viral clearance, the level of IL-6 rapidly decreases, weakening its inhibitory effect and further enhancing the net growth trend of immune suppression [S], which ultimately causes its steady state to remain slightly higher than the initial level.

From an immunophysiological perspective, the residual activation of humoral immunity [H] can be regarded as a “memory” response to potential reinvasion by pathogens, reflecting the persistence of humoral immunity. In contrast, the sustained high expression of immune suppression [S] may correspond to the regulatory system’s braking mechanism that suppresses excessive immune responses and maintains system stability. Thus, although the immune system tends to stabilize as a whole, its recovery process exhibits clear module asynchrony and functional delay. These results indicate that although viral clearance may be achieved within a short period after the termination of input, the full recovery of the immune system shows a pronounced delay effect, and the recovery speeds of different modules are not uniform. This dynamical feature is highly consistent with clinical observations of “residual immune activation during the recovery period after infection” (46).

To quantify the continuous durations of different phases during infection, two temporal indicators were defined over the observation interval [0, T_end], based on the state trajectories of [V] and IL-6 (47). As shown in Figure 6a, the thresholds θ_V and θ_IL−6 represent the critical values for infectiousness and host inflammatory response, respectively.

Figure 6

Infectious duration The onset time is defined as the earliest time when first reaches or exceeds the threshold . The offset time is defined as the earliest time when first falls below the threshold (). The infectious duration is then given by Equation 26.

Illness duration The onset time is defined as the earliest time when first reaches or exceeds the threshold . The offset time is defined as the earliest time when first falls below the threshold (). The illness duration is then given by Equation 27.

The numerical solutions were obtained using the built-in MATLAB solver ode15s, and the onset/offset times were determined from the original non-uniform time points by detecting the first crossing above or below the thresholds. To avoid spurious misjudgments caused by numerical noise, mild smoothing or a minimum duration window” can be applied, without altering the main form of the above definitions.

To further reveal how illness duration and infectious duration respond to key system parameters, Figure 6 illustrates the variations of these two dynamical indicators under multi-parameter perturbations. Figure 6a provides a schematic diagram of the infectious duration (gray shaded area) and the illness duration (blue shaded area), clearly defining the time intervals during which virus load [V] and IL-6 exceed their respective thresholds, and thereby allowing the calculation of infectious duration and illness duration. Figure 6 also presents the response characteristics of key system variables under parameter perturbations, as well as their sensitivity to the two physiological indicators, “illness duration” and “infectious duration.” The overall figure consists of seven subplots, each depicting different aspects of the immune system’s dynamical response under parameter perturbations. Figure 6a shows a typical virus–inflammation dynamical process. The upper panel illustrates the temporal evolution of the average viral concentration [V], while the lower panel shows the average level of the inflammatory factor IL-6 over time. The gray shaded region denotes the short-term infection phase, whereas the blue region denotes the phase of persistent infection or inflammatory activation. Around t ≈ 14days, the system undergoes viral input and reaches its peak, followed by transitions into different clearance or maintenance states. This figure is used to define the classification logic of the “illness period” and the “infectious phase”, virus load [V] peaks and subsequently falls below a certain threshold, while IL-6 remains at a relatively high level for a sustained period.

Figures 6b1–d2 present the variation curves of the key dynamical indicators, namely infectious duration and illness duration, with respect to different immune response–related parameters. In the calculations, each subplot examines the system’s response to the variation of a single parameter across different orders of magnitude, while all other parameters fluctuate randomly within ±30% probability. This design highlights both the independent effect of a specific parameter on system outputs and the intrinsic stochasticity of the system.

Specifically, to examine the effect of the viral self-replication half-saturation constant on infectious and illness durations, parameter scanning was performed on The results show that when is small (approximately less than 0.1), both infectious and illness durations are significantly shorter, indicating that viral replication can be effectively controlled by the immune system. However, once exceeds a critical threshold, the infectious and illness durations increase rapidly and tend toward a stable chronic infection state. This dynamical response highlights the pronounced nonlinear sensitivity of the system to viral replication parameters, where even small parameter changes can trigger dramatic shifts in system dynamics (Figures 6b1, b2).

In addition, we examined the effect of the parameter , which represents the efficiency of virus clearance mediated by humoral immunity. As increases, the ability of humoral immunity to clear the virus is enhanced, and both infectious duration and illness duration exhibit a gradually shortening trend. Unlike the sharp transitions induced by parameter , variations in lead to smoother and more continuous system responses. This result indicates that although humoral immunity–mediated viral clearance has a noticeable impact on system behavior, it is not a key parameter driving system switches or abrupt transitions (Figures 6c1, c2).

Furthermore, to analyze the effect of the parameter , which represents the strength of innate immunity self-activation, we compared its variation with the two parameters discussed above. Overall, changes in a_I1 exert a milder influence on infectious duration and illness duration. However, when is at relatively low levels (around 10⁻²), the system response still exhibits pronounced fluctuations and instability. This suggests that the self-feedback mechanism of innate immunity plays a strong nonlinear regulatory role in certain parameter regions, potentially inducing unstable infection states (Figures 6d1, d2).

In summary, the above analysis shows that parameters associated with different immune response mechanisms exert markedly different impacts on system dynamics. Among them, viral replication parameters (e.g., ) play a dominant role in controlling the switches of infection states, whereas variations in the humoral immunity parameter () and the innate immunity parameter () respectively reflect the secondary regulatory roles of viral clearance and immune modulation mechanisms in controlling infection duration. Such complex multi-module interactions of the immune system indicate that its overall function depends on both the coordination among modules and the timing and intensity of responses, thereby necessitating the consideration of dynamical sensitivity windows.

3.5 Different regulatory relationships give rise to distinct forms of system dynamics

Based on the above results, it can be seen that in addition to bistability, the system also exhibits pronounced oscillations and multistability under specific conditions (11, 48). Therefore, this section further focuses on the dynamical mechanisms underlying such complex behaviors. Specifically, we identify and analyze the core dynamical modules formed by the interactions among variables, and systematically investigate how key variable combinations influence these complex dynamical features. To this end, subsystems of equation (24) are studied separately, examining the interactions among variables within the system to gain deeper insight into the mechanisms by which oscillations arise.

First, considering that this study involves six key variables (viral load [V], innate immunity [I], cellular immunity [C], humoral immunity [H], immune suppression [S], and the inflammatory factor IL-6), we introduce subsystems for analysis. Specifically, each subsystem consists of a subset of these six variables, while the remaining variables not included are fixed at their steady-state levels. In this way, different combinations of the six variables yield a total of 2⁶ − 1 = 63 possible subsystem structures.

For each subsystem structure, numerical simulations were performed using single-parameter variation. Specifically, in all subsystems one key parameter was varied while all other parameters were kept fixed, and the dynamic behaviors were examined for each subsystem structure in turn. This approach enables a systematic investigation of whether different subsystem combinations can give rise to bistability or other complex dynamic behaviors.

As shown in Figure 7, we performed a frequency analysis of bistability occurrence in subsystems involving each variable. Group1 represents all subsystem configurations that include the given variable, while Group2 represents those that exclude it. Red indicates the frequency of bistability occurrence in subsystems involving the variable, whereas blue indicates the absence of bistability. The analysis reveals that different variables contribute to the emergence of complex system behaviors to varying degrees of importance. For example, subsystems involving virus ([V]), innate immunity ([I]), and cellular immunity ([C]) exhibit a notably higher proportion of bistability, suggesting that the interactions among these three variables are critical for the emergence of complex dynamic behaviors. In contrast, subsystems involving immune suppression ([S]), humoral immunity ([H]), and the inflammatory factor (IL-6) show relatively lower frequencies of bistability, indicating that interactions among these variables alone are insufficient to robustly generate complex dynamic features.

Figure 7

In summary, the above analysis identifies virus ([V]), innate immunity ([I]), and cellular immunity ([C]) as the core variable combination, which plays a critical role in shaping the system’s complex dynamic behaviors. This finding lays a solid foundation for further in-depth investigation of the dynamic mechanisms underlying virus–immune responses and provides a theoretical basis for the design of clinical intervention strategies.

4 Discussion

This study is based on a six-variable ordinary differential equation model incorporating viral load ([V]), innate immunity ([I]), cellular immunity ([C]), humoral immunity ([H]), immune suppression ([S]), and the inflammatory cytokine IL-6, to systematically investigate the nonlinear dynamic features of virus–immune system interactions. The model employs Hill functions to uniformly capture nonlinear activation, inhibition, and saturation effects across modules (22), and explicitly introduces an external input α(t) to characterize both continuous and short-term viral exposure scenarios. In this way, the immune network’s global response can be quantitatively assessed within a unified framework.

Bifurcation analysis of the model shows that, under continuous input, the relative magnitudes of viral replication rate and immune clearance efficiency directly determine the system’s steady-state landscape. When replication is strong and clearance is insufficient, the system may remain in a stable high-virus–high-inflammation state; when clearance is enhanced, the system can switch to a healthy low virus and low inflammation state. Two dimensional parameter scans further reveal that, under different parameter combinations, the system may exhibit monostability, coexistence of multiple steady states, or even strong sensitivity to initial conditions. This implies that even with similar exposure and parameter settings, long-term immune outcomes may differ substantially across individuals (41). Moreover, we identified the coexistence of “high/low states”, that is, even when viral levels are close to zero, cellular immunity and IL-6 can stabilize at distinct levels, indicating that positive and negative feedback can sustain multiple stable states under low viral pressure.

Under short-term input conditions, the system exhibits typicalmulti-timescale immune recovery dynamics, viral load declines rapidly once input ceases, while IL-6 simultaneously decays quickly, reflecting the transient nature of the inflammatory response. In contrast, humoral immunity and the immune suppression module remain elevated for a longer period and decrease more slowly, showing residual effects and asynchronous recovery. To quantitatively characterize this process, we defined the “infectious duration” and “illness duration” based on viral load and IL-6 thresholds, and examined their sensitivity under multiparameter perturbations. The results indicate that viral replication–related parameters play a dominant role in controlling the course of infection, capable of triggering nonlinear transitions from acute to chronic states, whereas humoral clearance efficiency and innate self-activation strength act as modulators, influencing the length and stability of infection and inflammation processes. These findings reveal that immune recovery after viral clearance is inherently unbalanced and delayed, providing a theoretical explanation for the clinically observed phenomenon of “residual immune activation during recovery”. Subsystem enumeration analysis further shows that the [V]–[I]–[C] triplet is the core structure driving multistability and complex oscillatory behavior, while [H] and the [S]–IL-6 module primarily affect the amplitude of responses and the synchronicity of recovery. Moreover, we identify the coexistence of distinct high and low immune activation states even when viral levels are close to zero, revealing that immune modules such as cellular immunity and IL-6 can stabilize at different levels due to internal feedback rather than direct viral drive. This finding provides a mechanistic basis for model reduction and the identification of potential intervention targets.

In summary, the mesoscale immune dynamics framework presented in this study elucidates how competition between viral replication and immune clearance generates multistability and sensitivity to initial conditions via bifurcation mechanisms, while emphasizing temporal disparities among immune modules during recovery. Temporal indicators are introduced to compare system outcomes across varying parameter settings and input modes. This framework enhances theoretical insights into virus-immune system coupling, furnishes methodological tools for parameter analysis and threshold identification, and yields practical guidance for immune intervention research. Future directions could involve fitting parameters with individualized longitudinal data, extending the model to incorporate stochasticity and spatial structures, and investigating regulatory strategies for key parameters to foster integration between mathematical modeling and clinical immunology.

References

• 1

  MurphyKWeaverC. Janeway’s immunobiology. New York: Garland science (2016).

  • Google Scholar

• 2

  GermainRN. Maintaining system homeostasis: the third law of newtonian immunology. Nat Immunol. (2012) 13:902–6. doi: 10.1038/ni.2404

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  DavisMMTatoCMFurmanD. Systems immunology: just getting started. Nat Immunol. (2017) 18:725–32. doi: 10.1038/ni.3768

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  MedzhitovRSchneiderDSSoaresMP. Disease tolerance as a defense strategy. Science. (2012) 335:936–41. doi: 10.1126/science.1214935

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  MalinziJJumaVOMadubuezeCEMwaonanjiJNkemGNMwakilamaEet al. Covid-19 transmission dynamics and the impact of vaccination: modelling, analysis and simulations. R Soc Open Sci. (2023) 10:221656. doi: 10.1098/rsos.221656

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  GoldblattD. SARS-CoV-2: from herd immunity to hybrid immunity. Nat Rev Immunol. (2022) 22:333–4. doi: 10.1038/s41577-022-00725-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  PepinKVolkovIBanavarJWilkeCGrenfellB. Phenotypic differences in viral immune escape explained by linking within-host dynamics to host-population immunity. J Theor Biol. (2010) 265:501–10. doi: 10.1016/j.jtbi.2010.05.036

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  PerelsonAS. Modelling viral and immune system dynamics. Nat Rev Immunol. (2002) 2:28–36. doi: 10.1038/nri700

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  ZitzmannCKaderaliL. Mathematical analysis of viral replication dynamics and antiviral treatment strategies: from basic models to age-based multi-scale modeling. Front Microbiol. (2018) 9:1546. doi: 10.3389/fmicb.2018.01546

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  GermainRNMeier-SchellersheimMNita-LazarAFraserID. Systems biology in immunology: a computational modeling perspective. Annu Rev Immunol. (2011) 29:527–85. doi: 10.1146/annurev-immunol-030409-101317

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  FerrellJ. Systems biology of cell signaling: recurring themes and quantitative models. New York: Garland Science (2021).

  • Google Scholar

• 12

  LiJWuJZhangJTangLMeiHHuYet al. A multicompartment mathematical model based on host immunity for dissecting covid-19 heterogeneity. Heliyon. (2022) 8:e09488. doi: 10.2139/ssrn.4021902

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  HunterCAJonesSA. Il-6 as a keystone cytokine in health and disease. Nat Immunol. (2015) 16:448–57. doi: 10.1038/ni.3153

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  TanakaTNarazakiMKishimotoT. Il-6 in inflammation, immunity, and disease. Cold Spring Harbor Perspect Biol. (2014) 6:a016295. doi: 10.1101/cshperspect.a016295

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  PawelekKAHuynhGTQuinlivanMCullinaneARongLPerelsonAS. Modeling within-host dynamics of influenza virus infection including immune responses. PloS Comput Biol. (2012) 8:e1002588. doi: 10.1371/journal.pcbi.1002588

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  PerelsonASNelsonPW. Mathematical analysis of hiv-1 dynamics in vivo. SIAM Rev. (1999) 41:3–44. doi: 10.1137/S0036144598335107

  • CrossRef
  • Google Scholar

• 17

  BaccamPBeaucheminCMackenCAHaydenFGPerelsonAS. Kinetics of influenza a virus infection in humans. J Virol. (2006) 80:7590–9. doi: 10.1128/JVI.01623-05

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  HanciogluBSwigonDClermontG. A dynamical model of human immune response to influenza a virus infection. J Theor Biol. (2007) 246:70–86. doi: 10.1016/j.jtbi.2006.12.015

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  SaenzRAQuinlivanMEltonDMacRaeSBlundenASMumfordJAet al. Dynamics of influenza virus infection and pathology. J Virol. (2010) 84:3974–83. doi: 10.1128/JVI.02078-09

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  BeaucheminCAHandelA. A review of mathematical models of influenza a infections within a host or cell culture: lessons learned and challenges ahead. BMC Public Health. (2011) 11:S7. doi: 10.1186/1471-2458-11-S1-S7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  CiupeSMHeffernanJM. In-host modeling. Infect Dis Model. (2017) 2:188–202. doi: 10.1016/j.idm.2017.04.002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22

  AlonU. An introduction to systems biology: design principles of biological circuits. New York: Chapman and Hall/CRC (2019).

  • Google Scholar

• 23

  WilkAJRustagiAZhaoNQRoqueJMartínez-ColónGJMcKechnieJLet al. A single-cell atlas of the peripheral immune response in patients with severe covid-19. Nat Med. (2020) 26:1070–6. doi: 10.1038/s41591-020-0944-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  HadjadjJYatimNBarnabeiLCorneauABoussierJSmithNet al. Impaired type I interferon activity and inflammatory responses in severe covid-19 patients. Science. (2020) 369:718–24. doi: 10.1126/science.abc6027

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25

  NowakMMayRM. Virus dynamics: mathematical principles of immunology and virology: mathematical principles of immunology and virology. UK: Oxford University Press (2000).

  • Google Scholar

• 26

  CaiYZhaoZZhugeC. The spatial dynamics of immune response upon virus infection through hybrid dynamical computational model. Front Immunol. (2023) 14:1257953. doi: 10.3389/fimmu.2023.1257953

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  SompayracLM. How the immune system works. Hoboken, NJ: John Wiley & Sons (2022).

  • Google Scholar

• 28

  AlbertsBJohnsonALewisJRaffMRobertsKWalterP. Molecular Biology of the Cell. 4 edn. New York: Garland Science (2002).

  • Google Scholar

• 29

  CainDWCidlowskiJA. Immune regulation by glucocorticoids. Nat Rev Immunol. (2017) 17:233–47. doi: 10.1038/nri.2017.1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  TrigunaiteADimoJJørgensenTN. Suppressive effects of androgens on the immune system. Cell Immunol. (2015) 294:87–94. doi: 10.1016/j.cellimm.2015.02.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  TangYLiuJZhangDXuZJiJWenC. Cytokine storm in covid-19: the current evidence and treatment strategies. Front Immunol. (2020) 11:1708. doi: 10.3389/fimmu.2020.01708

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  WherryEJBarouchDH. T cell immunity to covid-19 vaccines. Science. (2022) 377:821–2. doi: 10.1126/science.add2897

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  KaechSMCuiW. Transcriptional control of effector and memory cd8+ t cell differentiation. Nat Rev Immunol. (2012) 12:749–61. doi: 10.1038/nri3307

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  RastogiIJeonDMosemanJEMuralidharAPotluriHKMcNeelDG. Role of b cells as antigen presenting cells. Front Immunol. (2022) 13:954936. doi: 10.3389/fimmu.2022.954936

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  ZhouZLiDZhaoZShiSWuJLiJet al. Dynamical modelling of viral infection and cooperative immune protection in covid-19 patients. PloS Comput Biol. (2023) 19:e1011383. doi: 10.1371/journal.pcbi.1011383

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  CarsettiRZaffinaSPiano MortariETerreriSCorrenteFCapponiCet al. Different innate and adaptive immune responses to sars-cov-2 infection of asymptomatic, mild, and severe cases. Front Immunol. (2020) 11:610300. doi: 10.3389/fimmu.2020.610300

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  TakeuchiOAkiraS. Innate immunity to virus infection. Immunol Rev. (2009) 227:75–86. doi: 10.1111/j.1600-065X.2008.00737.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  KuznetsovYA. Elements of applied bifurcation theory. New York: Springer (1998).

  • Google Scholar

• 39

  DhoogeAGovaertsWKuznetsovYA. Matcont: a matlab package for numerical bifurcation analysis of odes. ACM Trans Math Software (TOMS). (2003) 29:141–64. doi: 10.1145/779359.779362

  • CrossRef
  • Google Scholar

• 40

  ChatterjeeBSingh SandhuHDixitNM. Modeling recapitulates the heterogeneous outcomes of sars-cov-2 infection and quantifies the differences in the innate immune and cd8 t-cell responses between patients experiencing mild and severe symptoms. PloS Pathog. (2022) 18:e1010630. doi: 10.1371/journal.ppat.1010630

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  EftimieRGillardJJCantrellDA. Mathematical models for immunology: current state of the art and future research directions. Bull Math Biol. (2016) 78:2091–134. doi: 10.1007/s11538-016-0214-9

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  MooreJBJuneCH. Cytokine release syndrome in severe covid-19. Science. (2020) 368:473–4. doi: 10.1126/science.abb8925

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  GuanW-jLiangW-hZhaoYLiangH-rChenZ-sLiY-met al. Comorbidity and its impact on 1590 patients with covid-19 in China: a nationwide analysis. Eur Respir J. (2020) 55:2000547. doi: 10.1183/13993003.00547-2020

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  KadelkaSBoribongBPLiLCiupeSM. Modeling the bistable dynamics of the innate immune system. Bull Math Biol. (2019) 81:256–76. doi: 10.1007/s11538-018-0527-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  XieQHuangDZhangSCaoJ. Analysis of a viral infection model with delayed immune response. Appl Math Model. (2010) 34:2388–95. doi: 10.1016/j.apm.2009.11.005

  • CrossRef
  • Google Scholar

• 46

  OpsteenSFilesJKFramTErdmannN. The role of immune activation and antigen persistence in acute and long covid. J Invest Med. (2023) 71:545–62. doi: 10.1177/10815589231158041

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  HadjichrysanthouCCauëtELawrenceEVegvariCDe WolfFAndersonRM. Understanding the within-host dynamics of influenza a virus: from theory to clinical implications. J R Soc Interface. (2016) 13:20160289. doi: 10.1098/rsif.2016.0289

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  NovákBTysonJJ. Design principles of biochemical oscillators. Nat Rev Mol Cell Biol. (2008) 9:981–91. doi: 10.1038/nrm2530

  • Pubmed Abstract
  • CrossRef
  • Google Scholar