A multicompartment mathematical model based on host immunity for dissecting COVID-19 heterogeneity.

The determinants underlying the heterogeneity of coronavirus disease 2019 (COVID-19) remain to be elucidated. To systemically analyze the immunopathogenesis of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection, we built a multicompartment mathematical model based on immunological principles and typical COVID-19-related characteristics. This model integrated the trafficking of immune cells and cytokines among the secondary lymphoid organs, peripheral blood and lungs. Our results suggested that early-stage lymphopenia was related to lymphocyte chemotaxis, while prolonged lymphopenia in critically ill patients was associated with myeloid-derived suppressor cells. Furthermore, our model predicted that insufficient SARS-CoV-2-specific naïve T/B cell pools and ineffective activation of antigen-presenting cells (APCs) would cause delayed immunity activation, resulting in elevated viral load, low immunoglobulin level, etc. Overall, we provided a comprehensive view of the dynamics of host immunity after SARS-CoV-2 infection that enabled us to understand COVID-19 heterogeneity from systemic perspective.

Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) causes coronavirus disease 2019 (COVID-19) and is responsible for the current pandemic. COVID-19 is recognized as a more complicated, multiorgan and heterogeneous illness than initially anticipated (Cordon-Cardo et al., 2020). Because of its heterogeneity, the spectrum of clinical features ranges from asymptomatic or mild upper respiratory tract symptoms to severe pneumonia and acute respiratory distress syndrome (Wu and McGoogan, 2020). COVID-19 patients display a diverse array of pathophysiological characteristics, including hyperinflammatory state, endothelial dysfunction and thromboembolic disease, as well as a clinical course that may be complicated by abrupt, unexpected deterioration during apparent recovery (Wiersinga et al., 2020). Several studies also found heterogeneity of immune response to vaccination in different individuals (Monin et al., 2021; Pimpinelli et al., 2021; Schwarz et al., 2021).

To investigate the heterogeneity of COVID-19, it is necessary to study immune response in different COVID-19 patients. SARS-CoV-2 infects cells expressing the surface receptor angiotensin-converting enzyme 2 (ACE-2) via the viral spike protein (Hoffmann et al., 2020). SARS-CoV-2 infection gives rise to host cell pyroptosis, damage-associated molecular pattern (DAMP) release and inflammatory cascade initiation (Vabret et al., 2020). A striking age-related disparity has been observed in the prevalence and severity of COVID-19, including differences in cross-neutralizing antibodies and differences in the levels and binding affinity of ACE-2 (Wong et al., 2020). The marked heterogeneity of disease prevalence and severity might be explained by age-dependent immunological mechanisms (Chen et al., 2021). Younger patients experience infrequent, mild, and self-limiting infections, possibly resulting from higher levels of cross-neutralizing antibodies, lower levels of ACE-2 receptors in nasal epithelium, plenty of naïve T/B cells and higher regulatory T-cell response, and lower IL-6 and TNF-α production from innate immune cells (Chen et al., 2021; Wong et al., 2020). Despite the many studies on heterogeneous immune response in COVID-19 patients, viral dynamics and host innate and adaptive immune responses in patients with diverse ages and health conditions, all in the setting of sparse antiviral and immunomodulating therapies, remain to be further studied.

Mathematical models have been utilized in the field of oncology to reconcile molecular reductionism with quantitative and holistic approaches that disentangle complex systems and provide a deeper understanding of cancer progression (Altrock et al., 2015; Anderson and Quaranta, 2008; Leon-Triana et al., 2021). Over the past two years, population-level epidemiological models helped predicting the spread of COVID-19 and provided insights for optimal control strategies (Wang et al., 2019; Zheng et al., 2021a, 2021b). Currently, mathematical modeling has also been used to analyze immune response in COVID-19 patients. Hernandez-Vargas et al. found a slow immune response against SARS-CoV-2 compared with influenza by fitting data (Hernandez-Vargas and Velasco-Hernandez, 2020). Jenner, A.L., et al. focused on innate immunity in different patients and proposed that delayed type-I IFN response caused severe tissue damage in the lungs, resulting in severe inflammation (Jenner et al., 2021). Voutouri, C., et al. incorporated both innate and adaptive immune responses, associating disease progression with the response rate of activated CD8+ T cells, and discussed the effect of multiple treatment strategies (Voutouri et al., 2021). Sadria and Layton simulated the actions of drugs that target SARS-CoV-2 virus infection and pointed out the importance of early intervention (Sadria and Layton, 2021). Although mathematical models have contributed in terms of host immunity against SARS-CoV-2 infection, we still need an immune model integrating multiple compartments in a manner that links peripheral blood and immunological processes, which is important for the detection of immune response (Ganusov and Auerbach, 2014). Thus, we constructed a multicompartment mathematical model to simulate the dynamics of host immunity in peripheral blood, lungs and secondary lymphatic organs (lymph nodes and spleen). We discussed the impact of key factors underlying the heterogeneity of COVID-19 patients, such as naïve T/B cell pools and the activation of antigen-presenting cells (APCs), in combination with the distributions of immune cells among the different compartments. Together, we provided a systematic framework for dissecting the heterogeneity of disease progression and symptomatology in patients with COVID-19.

Results

Construction of a simplified multicompartment mathematical model of immune response after SARS-CoV-2 infection

SARS-CoV-2 infection in an individual host will cause the activation of innate and adaptive immunity, including the proliferation of lymphocytes and the migration of immune cells and cytokines among lungs, draining lymph nodes, distant lymph nodes and spleen, and peripheral blood. In Figure 1, we illustrated the migration of immune cells and cytokines among the different compartments and the immune response network in lungs and secondary lymphatic organs.

Figure 1

Schematic illustration of multicompartment mathematical model and the host immune response network during the SARS-CoV-2 infection. (a) The model consisted of four compartments: draining lymph nodes, distant lymph nodes and spleen, peripheral blood and lungs. Chemokines/chemokine receptors regulate the migratory patterns and positioning of immune cells among the different compartments. (b) Compartment 3 (Lungs) was the primary site of immune response involving the interplay among SARS-CoV-2, bacteria, epithelial cells, innate cells and adaptive cells. Under the stimulation of both APCs and cytokines, SARS-CoV-2-specific naïve T/B cells proliferate and differentiate into effector T/B cells in secondary lymphoid organs (compartment 1 and 4), including cytotoxic T lymphocytes (CTLs), helper T (Th) cells, inducible regulatory T cells (iTregs) and antibody-secreting cells (ASCs). Th cells, including Th1, Th2, and Th17 cells, secrete various cytokines and facilitate immune cell function. nCoV (SARS-CoV-2 virus), Bac (bacteria), H (susceptible epithelial cells), I (infected epithelial cells), D (dead epithelial cells), Mono (monocyte), MΦ (macrophage), DC (dendritic cell), Neut (neutrophil), NK (natural killer cell). The adaptive immune cells pictured above are all virus-specific; non-virus-specific lymphocytes are not pictured.

During SARS-CoV-2 infection, SARS-CoV-2-infected lung epithelial cells first recruit innate immune cells, such as monocytes, NK cells, and neutrophils, to lungs via the secretion of chemokines. Monocytes are activated and differentiate into macrophages and dendritic cells (DCs), which work as antigen-presenting cells (APCs), upon exposure to viruses. Activated macrophages can phagocytize infected epithelial cells and virus, while activated DCs migrate to the secondary lymphoid organs to present antigens to activate naïve T cells. NK cells kill infected cells, and neutrophils phagocytize virus in the lungs.

Consequently, under the stimulation of both APCs and cytokines, naïve T/B cells proliferate and differentiate into SARS-CoV-2-specific effector T/B cells, including cytotoxic T lymphocytes (CTLs), helper T (Th) cells, inducible regulatory T cells (iTregs), and antibody-secreting cells (ASCs). CTLs migrate to the lungs and clear infected cells, and ASCs secrete immunoglobulin (Ig) which neutralizes virus. Th cells, including Th1, Th2, and Th17 cells, secrete various cytokines and facilitate immune cell function. iTregs exert immunosuppressive functions by secreting anti-inflammatory cytokines, including IL-10 and transforming growth factor β (TGF-β). MDSCs exhibit similar immunosuppressive functions, which can be observed in critically ill patients, especially critically ill patients with bacterial infection.

To investigate the clinical characteristics of patients with different severities, we analyzed longitudinal data of 194 COVID-19 patients from Wuhan Union Hospital, including their hemogram and serum cytokine profiles, as well as clinical classification and outcome. These clinical results suggested that the time courses of peripheral blood leukocytes, neutrophils, lymphocytes, IL-10 and IL-6 exhibited obvious distinguishability among mildly/moderately, severely, and critically ill patients (Figure S1). Particularly, the time courses of both peripheral blood lymphocytes and IL-6 are significantly related to the outcomes of COVID-19 patients (Figure S1). During the late stage of infection, critically ill patients showed significantly higher levels of IL-6, neutrophil count and MDSCs fraction than patients with mild and severe disease (Figure S1 and Figure S2), which could be attributed to elevated viral load and secondary bacterial infection in critically ill patients. Blood culture results confirmed bacterial co-infection in critically ill patients, who exhibited such common bacteria as klebsiella pneumoniae, staphylococcus epidermidis, and enterococcus faecium (Table S1).

Based on immunological principles, clinical data and COVID-19-related literature, we then constructed a multicompartment mathematical model to describe the dynamic processes of SARS-CoV-2 infection and the response of both innate immunity and adaptive immunity, together with leukocyte chemotaxis and post-viral bacterial infection. Details of the equations used in our model can be found in Section II of the Supplementary Material.

In our model, the dynamics of SARS-CoV-2-specific naïve T/B cells is depicted in the following equation, where [naïve] refers to the concentration of naïve T/B cells (109/L). In this article, naïve T/B cells mean SARS-CoV-2-specific naïve T/B cells without further specification.

We do not show intercompartmental flow and differentiation terms in the above equation. The kinetic parameters , and denote the proliferation rate, environmental bearing capacity and the apoptosis rate of naïve cells, respectively. describes the sufficiency of the naïve T/B cell pool, and larger indicates a more sufficient naïve cell pool.

In addition, as a major process during the initiation of immune response, antigen-bearing DCs and rare antigen-specific T/B cells must quickly find each other (Pulendran and Ahmed, 2011). Children displayed higher basal expression of relevant pattern recognition receptors, such as MDA5 (IFIH1) and RIG-I (DDX58), in upper airway epithelial cells. Consequently, macrophages, dendritic cells, and a pre-activated antiviral innate immunity could control SARS-CoV-2 infection early in children (Loske et al., 2021).

We utilize the following equation to describe APCs among which [Mono] and [I] denote the concentrations of monocytes and infected epithelial cells (109/L). Subscript 3 represents variation in the lungs (compartment 3).

We do not show intercompartmental flow terms and the effect of cytokines on activation in the above equation. Kinetic parameters and denote the differentiation rate and the apoptosis rate of APCs, respectively. We used the form of to describe the activation of APCs by SARS-CoV-2-infected epithelial cells, and K is inversely related to the sensitivity of APCs response to antigen.

For the sake of simplicity and clarity, we made the following main assumptions in our model.

First, we assumed that the parameters relevant to APCs function and naïve T/B cells decreased with disease severity. The characteristic of COVID-19 patients most commonly attributed to disease severity is age (Chen et al., 2021; Hu et al., 2021). This age dependency could be explained by the impaired immune response in elderly individuals with COVID-19, which was reported to be correlated with T/B cell repertoire restrictions and decreased DC antigen-presenting ability (Paschold et al., 2020; Zheng et al., 2020). Also, TCR/BCR repertoire composition and diversity have been considered as the major determinants of disease outcomes following viral infection (Gutierrez et al., 2020; Schultheiss et al., 2020).

Second, patients with severe COVID-19 have a higher viral load and a long virus-shedding period (Liu et al., 2020), which would cause even more severe damage to the lungs. Dysfunction of the mucosal immune system increases the susceptibility of patients to bacterial infection (Hanada et al., 2018), making severe patients more susceptible to bacterial infection than those who have not been seriously affected (Chen et al., 2020). In this regard, we integrated post-viral bacterial infection and clearance of bacteria by innate immunity into our model. We assumed that the growth rate of bacteria was linked to the reduced ratio of healthy lung epithelial cells, , where denotes the number of healthy lung epithelial cells, and denotes the initial value of .

Third, we assumed that high IL-6 level triggered the release of MDSCs. Severe patients with COVID-19 exhibit the emergence of myelopoiesis-generating immuno-suppressive myeloid cells (HLA-DR-/low monocytes and immature neutrophils) (Agrati et al., 2020; Giamarellos-Bourboulis et al., 2020; Silvin et al., 2020). MDSCs can be induced by inflammatory conditions, such as high concentrations of IL-6 (Giamarellos-Bourboulis et al., 2020; Schrijver et al., 2019; Tobin et al., 2019). Since host sepsis contributes to the extraordinary elevation of plasma IL-6 in critically ill patients with COVID-19 (Huang et al., 2020), factors generated during sepsis would induce the expansion and egress of MDSCs from the bone marrow into the peripheral blood (Schrijver et al., 2019).

Last, we assumed that high level of immunosuppressive cell-derived-TGF-β could trigger lymphocyte apoptosis (Banerjee et al., 2011). The concentration of serum TGF-β increased with time after SARS-CoV-2 infection (Ferreira-Gomes et al., 2021), which may suppress immune efficacy against SARS-CoV-2.

More detail about model assumptions and equations can be found in Section I and Section II of the Supplementary Material.

Modeling the dynamics of host immunity in different types of COVID-19 patients

Using the above model, we aimed to analyze and reveal the dynamic processes of virus and host immunity through replicating the migration of the lymphocytes among lungs, lymph nodes and peripheral blood in different patients after the onset of SARS-CoV-2 infection.

Under the constraints of the physiological level of immune cells and cytokines from clinical data and treatments, we estimated the ranges of kinetic parameters in the model to recapitulate the characteristics indicated by clinical data (See section Ⅳ in supplementary material for details). Accordingly, we classified COVID-19 patients into four typical disease severities: mild/moderate, severe, critically ill (survivor) and critically ill (non-survivor). To capture the clinical characteristics of patients with these four typical different severities, we selected and modified different sets of kinetic parameters, as shown in Table S2, especially different (environmental bearing capacity of naïve T/B cells) and (sensitivity of APCs response to antigen). We also decreased the killing rate of CTLs for severe and critically ill patients and lowered the clearance rate of bacteria for critically ill non-survivors. More detail about kinetic parameters can be found in Table S3. Initial state was set as steady state without viral infection (see Section III of the Supplementary Material).

First, in Figure 2, we longitudinally analyzed the dynamics of viral load, bacterial load, APCs, CTLs and SARS-CoV-2 specific Ig level in the lungs of COVID-19 patients. In the 1st week after SARS-CoV-2 infection, viral load in the lungs in the four types of patients showed no obvious difference. During the 2nd week after SARS-CoV-2 infection, in the lungs of the mild/moderate patients, APCs, CTLs and Ig level increased quickly. These activated innate/adaptive immune cells and Ig cooperated to clear the virus and kill infected epithelial cells, inducing the rapid decrease of viral load in mild/moderate patients (Figure 2a). However, in patients with higher severity, as shown in Figure 2b-d, during the 2nd week after SARS-CoV-2 infection, APCs, CTLs and Ig increased more slowly in comparison to the response in mild/moderate patients, in particular, for the critically ill survivor (Figure 2c) and non-survivor patients (Figure 2d). An obvious delayed activation of innate and adaptive immunity occurred and caused the accumulation of high viral load in the lungs and severe lung epithelial damage, setting up conditions for bacterial infection. In critically ill non-survivors during the 3rd and 4th week, bacterial infection did activate further innate immune response with a high level of APCs. Blood culture verified bacterial infections in some critically ill patients (Table S1). Our results indicated that the delayed response of APCs, CTLs and Ig contributed to a high viral load in the lungs and a long virus maintenance period in COVID-19 patients, both associated with higher severity.

Figure 2

Modelling the immune processes in lungs of four types of patients after SARS-CoV-2 infection. (a-d) Curves illustrate the dynamic processes of SARS-CoV-2 infection, bacterial infection, antigen-presenting cells (APCs), cytotoxic T lymphocytes (CTLs), and SARS-CoV-2-specific-immunoglobulin (Ig) level in the lungs (compartment 3).

Then, we investigated the dynamics of both IL-6 and lymphocytes, including all subsets of T cells and B cells, both SARS-CoV-2-specific and non-SARS-CoV-2-specific, and NK cells in the different compartments of four types of COVID-19 patients. We illustrated the time-dependent trajectories of IL-6 and lymphocytes in peripheral blood (Figure 3a-b), which qualitatively fit the clinical data for the four types of patients (Figure S1). Furthermore, in patients with higher severity, the simulation results not only showed a higher level of serum IL-6 (Figure 3a), but also suggested that extensive damage to lung epithelial, in turn, causes the accumulation of more macrophages (Figure S3). The simulation implied that different serum IL-6 level in COVID-19 patients could be explained by the difference in macrophage activation.

Figure 3

Simulation results of IL-6 and lymphocytes in four types of COVID-19 patients. (a-b) Dynamics of serum IL-6 and peripheral blood lymphocytes, including all subsets of T cells and B cells, both SARS-CoV-2-specific and non-SARS-CoV-2-specific, and NK cells. Different colored (line-style) lines represented different types of patients. Green, blue, purple, and red lines represented mild/moderate patients, severe patients, critically ill survivors and critically ill non-survivors, respectively. (c) Numbers of lymphocytes among the secondary lymphoid organs, peripheral blood and lungs in different types of COVID-19 patients.

Next, we investigated the mechanism of lymphopenia during SARS-CoV-2 infection. Recent studies show that the durability of lymphopenia determined disease severity (Tan et al., 2020). Disease severity is also correlated with the upregulation of chemokine and chemokine receptor expression (Bost et al., 2020; Chua et al., 2020), implying that severity-associated lymphopenia could be explained by the difference in lymphocyte chemotaxis. Moreover, the clinical study showed that abundance of inflammatory monocyte-derived macrophages increased in the lungs of critically ill COVID-19 patients (Liao et al., 2020; Yao et al., 2021).

In our model, we assumed that specific leukocyte trafficking molecules recruit leukocytes to SARS-CoV-2-infected lungs, while macrophages and infected epithelial cells are the main chemokine-secreting cells (Table S4). In Figure 3b, we showed that mild/moderate patients exhibited transient lymphopenia; severe and critically ill (survivors) patients had a longer period of lymphopenia; and critically ill (non-survivors) patients showed continuous lymphopenia. In Figure 3c, it can be seen that COVID-19 patients with different severities manifest the trafficking of lymphocytes among the secondary lymphoid organs, peripheral blood and lungs throughout the course of SARS-CoV-2 infection. In the early stage of SARS-CoV-2 infection (termed 0–10 days), we found that lymphocytes decreased in both peripheral blood and secondary lymphoid organs and were recruited to lung tissues and that more lymphocytes infiltrated into the lungs of COVID-19 patients with higher severity. In the recovered patients, following the resolution of inflammation, lymphocytes exited inflamed tissues and participated in recycling, together with the recovery of the total number of lymphocytes. However, critically ill patients (non-survivors) at the late stage of disease (termed 21–28 days) exhibited the continuous decrease of lymphocytes in peripheral blood and secondary lymphoid organs, but a high level of lymphocytes in the lungs (Figure 3c).

We also presented the time courses of immune cells and cytokines in the different compartments, as well as viral load, bacterial load, and epithelial cells in the lungs in Figures S3, S4, S5, S6 and S7 where we illustrated the disparate immune response among the different types of patients in the different compartments.

Furthermore, we performed parametric sensitivity analysis of our model and listed the top 25 (≈11%) parameters according to their sensitivities in Figure S8, where , , parameters relevant to viral infectivity (e.g., the infection rate of SARS-CoV-2 , apoptosis rate of infected epithelial cell ), etc., showed high sensitivity to the change of values of four typical disease severities (Table S2). In this paper, we focused on the heterogeneity of immune system in different patients, and did not take the infectivity of different viurses into account. Other parameters with high sensitivity, such as apoptosis rate of neutrophils , the volume of compartment 4 , were not likely to vary much in different patients, so we did not modify them. We also simulated individual difference within the same patient type by perturbation of all parameters, as shown in Figure S9. Although individual differences of immune response within the same patient type exist, we still focused on the heterogeneity among the four typical patient types in this article.

Activation and strength of both innate and adaptive immunity determined COVID-19 severity

Based on our recent work (Zhou et al., 2021), we developed quantitative indicator to describe the ability and strength of host innate and adaptive immunity against SARS-CoV-2 infection. We analyzed the two equations related to SARS-CoV-2 virus and infected epithelial cells. In the lungs of patients, , and denote the concentrations of virus, healthy epithelial cells and infected epithelial cells, respectively. Thus, we have:Where, , , , , , and denote the concentration of Type-Ⅰ interferon, neutrophil, immunoglobulin, macrophage, NK, and CTL respectively. The detailed information of parameters is listed in Table S3. Assuming that the infected epithelial cells rapidly reach equilibrium, , we simplified Eq. (3) to

We defined as the immune strength of host that consists of both innate and adaptive immunity. Similarly, we also respectively defined and as innate and adaptive immune strength,

We plotted the time courses of and in patient lungs, as shown in Figure 4a. In mild/moderate patients, we see the quickest activation of innate immune response and highest level of adaptive immune strength. For patients with more severe disease, the initial activation time of innate and adaptive immunity occurred later. A longitudinal analysis of , and can be found in Figure S10.

Figure 4

Connecting lymphocytes and IL-6 levels in peripheral blood with the strength of patients' immune response. (a) Simulation results of innate immune strength and adaptive immune strength in COVID-19 patients with different disease severities. (b) The dynamic trajectories of patients with different severities are plotted on the innate immune strength and serum IL-6 plane. (c) The dynamic trajectories of patients with different severities are plotted on the adaptive immune strength and peripheral blood (PB) lymphocytes plane. The (b) and the (c) were normalized to the maximum value for better comparison. The yellow circle (b and c) indicated the starting point of SARS-CoV-2 infection. The arrows on the curve (b and c) indicated the direction of the trajectory, and the red highlight indicated bacterial infection. (d) The level of serum IL-6 changed with the different levels of bacterial load. Different levels of bacterial load were simulated by modifying the parameter (infection rate of bacteria) based on the parameters of critically ill survivors, where was set as 0, 5, and 80 (1/day) to represent no, low and high bacterial loads, respectively. (e) The total number of lymphocytes in the secondary lymphoid organs, peripheral blood and lungs changed with different levels of MDSCs. We set the parameter (Hill coefficient of IL-6 enhancing the release of MDSCs) to 1000, 400, and 100 (pg/ml) to simulate low, median and high MDSCs levels, respectively.

Next, we connected host immune response status to the levels of serum IL-6 and peripheral blood lymphocytes. We found that serum IL-6 was positively associated with innate immune strength , which then means that cytokine release syndrome correlated with the strong activation of innate immunity (Figure 4b). In critically ill patients, secondary bacterial infection also caused strong activation of innate immunity, leading to a high level of serum IL-6 during the 3rd and 4th weeks. Based on the critically ill survivor model, we changed bacterial load in the lungs by setting to 0, 5, and 80 (1/day) to represent low, median and high bacterial loads in the lungs, respectively. In Figure 4d, we found that the level of serum IL-6 was significantly elevated when bacterial load increased in the lungs. In particular, a low bacterial load could cause a second peak of serum IL-6 level, suggesting that bacterial infection can further exacerbate cytokine storm.

In addition, the recovery of peripheral blood lymphocyte level served as a positive sign of establishment of adaptive immunity, which is indicated by . In the mild/moderate and severe patients (upper row of Figure 3c and Figure 4c), when began to increase, the establishment of adaptive immunity could protect lung tissues from viral damage and decrease the infiltration of lymphocytes into the lungs, thereby facilitating the recovery of lymphocytes in both peripheral blood and secondary lymphoid organs. However, from the lower row of Figure 3c and Figure 4c, we found that low level of and secondary bacterial infection led to the sequestration of lymphocytes in the lungs in critically ill patients, inhibiting the recovery of lymphocytes in peripheral blood.

In Figure 4e, we simulated the effect of different MDSCs levels on the level of lymphocytes by setting the Hill coefficient of IL-6 to enhance the release of MDSCs to 1000, 400, and 100 (pg/ml), and we found MDSCs to be negatively correlated with the total number of lymphocytes in all compartments. TGF-β derived by MDSCs may induce lymphocyte apoptosis and cause the continuous lymphopenia in critically ill patients. Our simulation showed that proliferation and natural apoptosis in different patient types were similar; however, TGF-β-induced apoptosis was higher in critically ill patients than that in other patients (Figure S11). Thus, we revisited Figure 3b-c and found that the total number of lymphocytes significantly decreased in critically ill non-survivors, which included a decrease in lymphocyte numbers in peripheral blood and secondary lymphoid organs. By contrast, the total number of lymphocytes almost recovered at the 28th day after infection in critically ill survivors. Given that abundant MDSCs were enriched in myeloid cells in critically ill patients (Figure S2), we concluded that continuous lymphopenia in critically ill patients might result from MDSCs inhibiting lymphocyte activity and further inducing lymphocyte apoptosis.

Naïve T/B cell pools and activation of APCs impact the clinical outcome of patients

To reveal the roles of naïve T/B cell pools and the activation of APCs, we theoretically investigated how the change of and influences the dynamic processes of immune response and clinical outcome in COVID-19. For critically ill (survivors) patients, we found in Figure 5a-c that increasing levels of both and contributed to the decrease of average viral load, bacterial load and serum Ig in the lungs, suggesting a better clinical outcome. The convexity of contours also implied a synergistic role for and . We further illustrated the impact of and in Figure S12, respectively. High played a role in establishing a strong adaptive immunity at a relatively low innate immunity; high played a role in the control of viral growth in the early phase and triggered the quick establishment of strong adaptive immunity (Figure S12).

Figure 5

SARS-CoV-2-specific naïve T/B cell pools and activation of APCs impact the clinical status and outcome of patients. (a-c) represents the environmental carrying capacity of antigen-specific T/B naïve cells, and reflects the activation threshold of APCs. Here, and were multiples relative to the parameters of the critically ill (survivor) type denoted by purple star. Analysis detailing how and affected the following variables: (a) average viral load in the lungs, (b) average bacterial load in the lungs, and (c) average immunoglobulin (Ig) level in peripheral blood (PB). The red line in (b) indicates the boundary of bacterial infection. Based on the corresponding parameter values = 0.7 = 0.7 (red square at bottom left in (b)), = 1 = 1 (purple star at middle in (b)), and = 1.5 = 1.5 (blue triangle at upper right in (b)), we calculated and plotted the pseudopotential and the typical trajectories of lymphocytes and IL-6 (d-f). Red dashed lines show the contour of a specific energy for better visualization. For high levels of and , patients exhibited low levels of viral load, bacterial load and serum Ig level; for low levels of and , patients exhibited severe bacterial infection and high levels of MDSCs (state A). By increasing the fold change of and , the barrier obstructing the system from recovery disappeared, and state A attractor also disappeared as no bacterial infection occurred.

Furthermore, we changed and based on the parameters of critically ill patients (survivor) (Table S8) and found the pseudopotential perspective of three typical COVID-19 trajectories (Figure 5d-f). In Figure 5d, when the system had relatively low values of and , it would start from the initial state and evolve to stable state A, or low lymphocyte level and high IL-6 level in peripheral blood, that corresponded to critically ill (non-survivor) patients with cytokine storm. These patients typically had severe bacterial infections at the late stage, as shown in the lower left corner of Figure 5a-c. When and increased, as shown in Figure 5e and Figure 5f, state A became unstable. In Figure 5f, bacterial infection did not occur, and the system evolved back to the healthy state with a high lymphocyte level and low IL-6 level in peripheral blood. Similarly, we illustrated in Figure S13 the pseudopotential landscape of the mild/moderate, severe, critically ill (survivor), and critically ill (non-survivor) patients.

Thus, our model quantitatively illustrated that the naïve T/B cell pools and the activation of APCs had a significant impact on the clinical outcome of COVID-19 patients. During the interaction between SARS-CoV-2 and host immunity, the hosts with quick activation of innate and adaptive immunity and sufficient naïve T/B cell pools would quickly recover to health; however, the hosts with diminished innate and adaptive immunity activation and insufficient naïve T/B cell pools would most likely have severe disease and bad clinical outcome.

Discussion and conclusion

In this study, we constructed a multicompartment mathematical model to assess underlying immunological mechanisms, clinical progression and disease severity, aiming to analyze the heterogeneity of COVID-19 patients. Based on the model, we simulated and predicted the dynamics of host immunity after SARS-CoV-2 infection and obtained the following results. (1) Early-stage lymphopenia was related to lymphocyte chemotaxis. (2) The prolonged lymphopenia in critically ill patients was associated with MDSCs, which were induced by highly inflammatory environment especially during the bacterial infections. A recent study also indicated that the numbers of suppressive immature neutrophils and/or G-MDSCs expanded during severe COVID-19 infection, which was associated with lymphopenia and disease severity (Penaloza et al., 2021). (3) The delayed onset of innate and adaptive immune responses after SARS-CoV-2 infection would cause a series of cascading reactions, including elevated viral load in the lungs and secondary bacterial infection, in the absence of medical interventions. Secondary bacterial infection will subsequently evoke a strong innate immunity, leading to cytokine storm and high levels of MDSCs. (4) SARS-CoV-2-specific naïve T/B cell pools and activation of APCs had a significant impact on COVID-19 severity. Patients with insufficient SARS-CoV-2-specific naïve T/B cells and inactive APCs would likely lapse into continuous lymphopenia and cytokine storm.

Our results, especially the simulation results in Figure 5, suggest that the individual with rapid and strong activation of immunity will quickly recover to health. Rapid and strong immune activation relies on active APCs and sufficient SARS-CoV-2-specific naïve T/B cell pools. This provide some hints for the further clinical treatments. On the one hand, enhancing APCs activation during the early phase of infection may be a treatment option for high-risk COVID-19 patients. Recent study proposed that early intervention with recombinant interferon-α2b helps reduce both viral replication and secondary viral infection of neighboring cells, thus providing a rational treatment regimen for the management of COVID-19 (Pandit et al., 2021). On the other hand, vaccination induces the SARS-CoV-2-specific memory T cells and memory B cells. Memory T/B cells will rapidly proliferate and differentiate into effector T/B cells after SARS-CoV-2 infection, which provides a quick and strong response of adaptive immunity. This will increase efficacy against severe disease and reduce mortality. Notably, the T/B cell pools here refer to the virus-specific T/B cell pools, and variants may evade vaccine-induced response. Thus, it is necessary to enlarge memory T/B cells in terms of quantity and coverage, and combining different COVID-19 vaccines may be a helpful strategy against variants. In our recent unpublished work, we developed a mathematical model including memory T/B cells and further discussed the efficacy of vaccines against different variants.

In summary, our mathematical model provides a framework to understand the dynamics of immune processes, enabling us to reveal the immunopathogenesis of COVID-19. Further, this mathematical framework may be useful for exploring the effects of medication and vaccines on the immune process as well as optimizing treatment strategies accordingly. It may also be applied to study immunodynamics of other diseases.

Limitations of the study

Our mathematical model has some limitations. First, although the results of our model qualitatively matched our clinical observations, the lack of detailed clinical data in both the lungs and secondary lymphatic organs makes it unable to fit parameters in our model. Therefore, more clinical data, especially animal model data, are needed for precise simulation and verification. Second, we mainly focused on how naïve T/B cell pools and APC activation affect severity, yet other factors may also affect severity, such as immune cell function and immunosuppression (see Section Ⅳ in Supplementary Material). Third, to simplify our model, adaptive immunity triggered by bacterial infection or other pathogenic infections was not included in our model. Based on our mathematical model, future model development can include additional complexities of host immune response against variant strains, as well as modulation of the immune system by vaccination.

Materials and methods

Basic framework of the mathematical model

The model contains four compartments: draining lymph nodes, peripheral blood, lung and distant lymph nodes and spleen. The model consists of epithelial cells, SARS-CoV-2, bacteria, cytokines, innate immune cells (neutrophils, monocytes, macrophages, dendritic cells (DCs), myeloid-derived suppressor cells (MDSCs) and natural killer (NK) cells), and adaptive immune cells (T cell and B cell subsets). Chemokines/chemokine receptors control the migratory patterns and positioning of immune cells among the different compartments. We used ordinary differential equations (ODEs) to simulate the time-dependent functions of immunologic variations in the different compartments. Our model consists of equations for 109 immunologic variations (Table S5), which contain 223 parameters (Table S3). According to the migration patterns of immune cells and cytokines, we categorized them into six groups (Table S6). Colony-stimulating factors (CSFs) promoted the mobilization of monocytes and neutrophils from the bone marrow into the peripheral blood (Table S7). In this work, we used Hill-type functions to quantitatively describe the activation of APCs and the interactions between immune cells and cytokines. ODEs were encoded with Python 3.7, and the odeint from the scipy.integrate package was used to solve the ODEs. Detailed information on equations and parameters is included in the Supplementary Material.

Model validation with clinical data of various COVID-19 phenotypes

A total of 194 laboratory-confirmed COVID-19 admitted cases with clarified outcome, either discharged or deceased, were collected at Union Hospital of Tongji Medical College in Huazhong University of Science and Technology (Wuhan, P. R. China) (Table S1). Clinical information for all patients was collected from the hospital's electronic history system. Of the 194 patients, 42 were utilized for the detection of immune cell subsets by flow cytometry during the period of lymphopenia (11 mild/moderate cases, 12 severe cases and 19 critically ill cases). Blood samples collected at any time from mildly/moderately ill patients were used as control groups because these patients experienced transient lymphopenia or no lymphopenia. Disease severity was assessed according to the Seventh Version of the Novel Coronavirus Pneumonia Diagnosis and Treatment Guidance from the National Health Commission of China (Diagnosis and Treatment Protocol for Novel Coronavirus Pneumonia (Trial Version 7), 2020). This study was conducted in accordance with the Declaration of Helsinki and was approved by the Ethics Committee of Union Hospital, Tongji Medical College, Huazhong University of Science and Technology (#2020/0004). Written informed consent was waived owing to the emergence of this high-risk infectious disease.

Method for analyzing the impacts of the naïve T/B cell pools and APCs on clinical status

is a set of parameters, including , which represent the environmental carrying capacity of SARS-CoV-2-specific naïve CD4+ T, naïve CD8+ T, and naïve B cells in compartment 1 and compartment 4. represents the Hill coefficient of APCs differentiation. Multiples of multiply , and the multiples of divide . The average viral load, the average bacterial load, and the average Ig level were defined as:

Here, represents the end time of the simulation, which is the 28th day in our model. The peak time of viral load refers to the time corresponding to the peak viral load.

Method of pseudo landscape

We referred to the method of Zhou J X, et al. (Zhou et al., 2012) to calculate the pseudo landscape to better visualize the effect of APC capacity and the SARS-CoV-2-specific naïve T/B cell pools on clinical conditions. The pseudo potential is defined as:where represents the probability that the trajectory passes through , which is defined as ; is the number of times that the trajectory passes through and is a technical approximation to prevent divergent results. We divided the plane of into a 128∗128 grid. To obtain the different trajectories of COVID-19, we added the disturbances to every parameter within the range of [90%,110%] using the Latin hypercube sampling scheme. After each step of the simulation ( = 0.25 days), the output is multiplied by Gaussian noise (μ = 1, σ = 1/20) and then used as the initial value of the next step. We computed 16384 different trajectories and excluded those with problems (e.g., nan appears in the output trajectory from random noise); we then performed statistics on the trajectories. Gaussian blur (windows of 3∗3) was used for noise reduction.

Declarations

Author contribution statement

Jianwei Li, Jingpeng Zhang: Analyzed and interpreted the data; Wrote the paper.

Jianghua Wu: Conceived and designed the experiments; Performed the experiments; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Lu Tang: Performed the experiments.

Heng Mei, Yu Hu: Conceived and designed the experiments.

Fangting Li: Conceived and designed the experiments; Wrote the paper.

Funding statement

This work was supported by the National Key R&D Program in China (Grants No. 2020YFC0845700, No. 2018YFA0900200 and No. 2020YFA0906900), the National Natural Science Foundation of China (Grants No. 12174007 and No. 12090054), Scientific Research Projects of the Chinese Academy of Engineering (Grant No. 2020-KYGG-01-07) and Fundamental Research Funds for the Central Universities (Grant No. 2020kfyXGYJ029).

Data availability statement

Data included in article/supplementary material/referenced in article.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.