https://orcid.org/0000-0001-9103-7092
Figures
Abstract
To understand the diversity of immune responses to SARS-CoV-2 and distinguish features that predispose individuals to severe COVID-19, we developed a mechanistic, within-host mathematical model and virtual patient cohort. Our results suggest that virtual patients with low production rates of infected cell derived IFN subsequently experienced highly inflammatory disease phenotypes, compared to those with early and robust IFN responses. In these in silico patients, the maximum concentration of IL-6 was also a major predictor of CD8+ T cell depletion. Our analyses predicted that individuals with severe COVID-19 also have accelerated monocyte-to-macrophage differentiation mediated by increased IL-6 and reduced type I IFN signalling. Together, these findings suggest biomarkers driving the development of severe COVID-19 and support early interventions aimed at reducing inflammation.
Author summary
Understanding of the diversity of immune responses to SARS-CoV-2 infections is critical for improving diagnostic and treatment approaches. Identifying which immune mechanisms lead to divergent outcomes can be clinically difficult, and experimental models and longitudinal data are only beginning to emerge. In response, we developed a mechanistic, mathematical and computational model of the immunopathology of COVID-19 calibrated to and independently validated against a broad set of experimental and clinical immunological data. To study the drivers of severe COVID-19, we used our model to expand a cohort of virtual patients, each with realistic disease dynamics. Our results suggest key processes that regulate the immune response to SARS-CoV-2 infection in virtual patients and suggest viable therapeutic targets, underlining the importance of a rational, multifaceted approach to studying novel pathogens using intra-host models.
Citation: Jenner AL, Aogo RA, Alfonso S, Crowe V, Deng X, Smith AP, et al. (2021) COVID-19 virtual patient cohort suggests immune mechanisms driving disease outcomes. PLoS Pathog 17(7): e1009753. https://doi.org/10.1371/journal.ppat.1009753
Editor: Andrew J. Yates, Columbia University Medical Center, UNITED STATES
Received: April 22, 2021; Accepted: June 24, 2021; Published: July 14, 2021
Copyright: © 2021 Jenner et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the manuscript and its Supporting information files.
Funding: ALJ was supported by Fonds de recherche du Québec-Santé Programme de bourse de formation postdoctorale pour les citoyens d’autres pays and a Centre for Applied Mathematics in Bioscience and Medicine (CAMBAM) Postdoctoral fellowship: http://www.frqs.gouv.qc.ca, http://www.crm.umontreal.ca/labo/cambam/en/. VC was supported by an IVADO Undergraduate Introduction to Research Scholarship. RA and AMS were supported by NIH R01 AI139088: https://www.nih.gov. ALJ and MC were supported by NSERC Discovery Grant RGPIN-2018-04546 and NSERC Alliance COVID-19 ALLRP 554923 – 20: https://www.nserc-crsng.gc.ca/index_eng.asp. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Clinical manifestations of SARS-CoV-2 infection are heterogeneous, with a significant proportion of people experiencing asymptomatic or mild infections that do not require hospitalization. In severe cases, patients develop coronavirus disease (COVID-19) that may progress to acute respiratory distress syndrome (ARDS), which is frequently accompanied by myriad inflammatory indicators [1]. Mounting evidence points to a hyper-reactive and dysregulated inflammatory response characterized by overexpression of pro-inflammatory cytokines (cytokine storm) and severe immunopathology as specific presentations in severe COVID-19 [2–6]. An over-exuberant innate immune response with larger numbers of infiltrating neutrophils [7,8] arrests the adaptive immune response through the excessive release of reactive oxygen species that leads to extensive tissue damage and depletion of epithelial cells [9]. In addition, T cell lymphopenia, in particular, is one of the most prominent markers of COVID-19 and has been observed in over 80% of patients with severe disease [6,10–12]. However, the immune mechanisms that lead to disparate outcomes during SARS-CoV-2 infection remain to be delineated.
Cytokines are critically important for controlling virus infections [13,14] and are central to the pathophysiology of COVID-19, sometimes playing a detrimental role in the context of a cytokine storm [10]. For example, interleukin-6 (IL-6) can stimulate CD8+ T cell expansion under inflammatory conditions [15]; however, in hospitalized SARS-CoV-2 patients with lymphopenia, IL-6 has been shown to be elevated [16] without an increase in CD8+ T cell counts [17]. Type I interferons (such as IFNs- α, β [18]) also play a major role in limiting viral replication by inducing a refractory state in susceptible and infected cells [19–21]. Due to this, it has been suggested that a delay in mounting an effective IFN response may be responsible for COVID-19 severity [10] as it is for other highly pathogenic coronavirus (i.e. SARS-CoV and MERS) infections [13]. Interferon delays/dysregulation contributing to COVID-19 severity are a known feature of coronavirus infections in humans [22] and type-I interferons are known to downregulate viral replication in infected and neighbouring cells [22]. The emerging experimental and clinical pictures of the impact of IFN-I dynamics on SARS-CoV-2 infections is consistent with both SARS-CoV and MERS [23–26]. Further, retrospective and clinical trials have shown that earlier administrations of IFN was associated with or significantly improved survival rates [27,28]. Overall, patients with severe COVID-19 present with lymphopenia [14,29], and are likely to have increased inflammatory cytokines such as IL-6, granulocyte-macrophage colony-stimulating factor (GM-CSF), and granulocyte colony-stimulating factor (G-CSF) [7,17,30].
Because the identification of immune mechanisms responsible for divergent disease outcomes can be difficult clinically, and experimental models and longitudinal data are only beginning to emerge, theoretical explorations are ideal [31]. Quantitative approaches combining mechanistic disease modelling and computational strategies are being increasingly leveraged to investigate inter- and intra-patient variability by, for example, developing virtual clinical trials [32–34]. Such in silico trials enable the theoretical exploration of how multiple key underlying mechanisms simultaneously impact disease severity [34]. More recently, viral dynamics models [35,36] have been applied to understand SARS-CoV-2 within-host dynamics and their implications for therapy [37–42]. However, there are few comprehensive models that integrate detailed immune mechanisms and allow interrogation of the dynamics controlling divergent outcomes, and none have attempted to quantify the high degree of variability in patient responses to SARS-CoV-2 through modelling.
In this study, we developed a mechanistic mathematical model to describe the within-host immune response to SARS-CoV-2. We explicitly modelled the interactions between epithelial cells, innate and adaptive immune cells, and cytokines. The model was fit to various in vitro, in vivo, and clinical data, analyzed to predict how early infection kinetics facilitate downstream disease dynamics, and used to create a virtual patient cohort with realistic disease courses. Our results suggest that mild and severe disease are distinguished by the rates of monocyte differentiation into macrophages and of IFN production by infected cells. In our virtual cohort, we found that severe COVID-19 responses were tightly correlated with a delay in the peak IFN concentration and that a large increase in IL-6 was the dominant predictor of CD8+ T cell depletion. These results provide insight into differential presentations of COVID-19 by suggesting key regulators of severe disease manifestation particularly related to monocyte differentiation and IL-6 concentrations. Importantly, these key mechanisms suggest promising avenues of experimental research into the mechanisms of immunopathology in COVID-19.
Results
Modelling the immune response to SARS-CoV-2 and the impact of delayed IFN on infection dynamics
To study the dynamics of SARS-CoV-2 infection and the development of COVID-19, we constructed a computational biology model of host-pathogen interactions (Eqs. S1-S22, with variables and parameters summarized in S1 Table and schematic in Fig 1). The model includes susceptible lung epithelial cells (S) that encounter virus (V) and become infected (I) before turning into damaged or dead cells (D) due to viral infection or immune involvement. The immune response is orchestrated by cytokines that act to stimulate the immune cell subsets present in the tissues and recruit cells from the bone marrow and circulation (Fig 1A). We have chosen to model type I IFN, IL-6, G-CSF, and GM-CSF because of their roles in viral resistance, inflammation, recruitment and differentiation, respectively. Upon infection, lung epithelial cells secrete type I IFNs (F) that cause adjacent cells to become resistant to infection (R) and decrease the production of newly infected cells [43]. Through stimulation by infected and dead cells, alveolar (lung tissue-resident) macrophages (MΦR) become inflammatory macrophages (MΦI), which also arise through monocyte (M) differentiation following stimulation by GM-CSF (G) or IL-6 (L) [44]. Neutrophils (N) are recruited to the infection site by G-CSF and release reactive oxygen species (ROS) causing bystander damage to infected and susceptible cells [45,46]. CD8+ T cells (T) are subsequently recruited to the infection site following a delay to account for antigen presentation, with expansion modulated by type I IFN and IL-6 concentrations. See Materials and methods for a complete description.
The model in Eqs. S1-S22 reduced to A) cell dynamics B) cytokine production dynamics and C) cytokine binding kinetics. Unique lines represent induced cell death (double line), recruitment (dashed line), cell type change or production (solid line), and cytokine production (square arrow). Cell and/or cytokines along joining lines denote a causal interaction. A) Virus (V) infects susceptible lung epithelial cells and creates either infected (I) or resistant (R) cells depending on the concentration of type I IFN. Infected cells then either die and produce new virus or are removed via inflammatory macrophages (MΦI) or CD8+ T cells (T) that induce apoptosis to create dead cells (D). Neutrophils (N) cause bystander damage (death) in all epithelial cells and are recruited by individually G-CSF and IL-6 concentrations. CD8+ T cells are recruited by infected cells and their population expands from IFN signalling. T cell recruitment is inhibited by IL-6 concentrations. Monocytes (M) are recruited by infected cells and GM-CSF and differentiate into inflammatory macrophages based on the individual concentrations of GM-CSF and IL-6. Tissue-resident macrophages (MΦR) also become inflammatory macrophages through interaction with dead and infected cells. Dead cells are cleared up by inflammatory macrophages and also cause their death. B) Type I IFN is produced by infected cells, inflammatory macrophages and monocytes. G-CSF is produced solely by monocytes and GM-CSF is produced by monocytes and macrophages. IL-6 is produced by monocytes, inflammatory macrophages and infected cells. C) Cytokine receptor binding, internalization and unbinding kinetics considered for each cell-cytokine interaction.
Because the model has several parameters that are undetermined biologically and insufficient data exists to confidently estimate their values, we used a stepwise approach to parameter estimation (see Materials and methods and S1–S5 Figs). We chose this approach to ensure that individual aspects of a given biological interaction could be better understood in isolation and to help reduce the degrees of freedom at each stage of parameter fitting. We first confirmed that we could recapitulate early infection viral kinetics with a reduced version of the full model (‘viral model’). For this, we excluded immunological variables (i.e. only including Eqs 6–9) and estimated parameters relating to viral kinetics by fitting to viral load data from macaques [47] and then viral load in humans hospitalized in Singapore [48] and Germany [49] (see Materials and methods). The resulting model dynamics were in good agreement to these early infection data (see S6 Fig for the macaque and Fig 2 for the human SARS-CoV-2 viral dynamics) and demonstrate a rebound in epithelial lung tissue as the viral load and infected cells decrease.
A reduced version of the full model (all cytokine and immune cells set to 0, Eqs 6–9) was fit to data from hospitalized patients, after initial estimation from viral loads in macaques [47] (S6 Fig) to estimate preliminary viral kinetic parameters. A) Virus (V) infects susceptible cells (S) making infected epithelial cells (I) which then die to produce dead cells (D) and new virus. B) Viral load data (log10(copies/mL) from eight human patients (three from Singapore S5, S6 and S18, and five from Germany, G1, G2, G5, G6, G7) were digitized from previous results [50], and parameters from the viral dynamics submodel were estimated using a non-linear squares optimization routine. β, dI, V0 and dV were estimated from the reduced viral dynamics model in A) (see Methods and S1 Table). Individual patient measurements are depicted by coloured circles. Solid black line: average model prediction; grey shaded region: predicted standard deviation from average. S (time axis) indicates the day of symptom onset.
We then isolated the IFN dynamics to assess clinical and experimental findings suggesting that delaying IFN results in more severe presentations in highly pathogenic coronavirus infections including SARS-CoV-2 [10,13,14]. Using the parameters obtained from the ‘viral model’ (Eqs 6–9; S1 Table), we then simulated the impact of IFN with the ‘IFN model’ (Eqs 10–16 and Fig 3). We examined the predicted dynamics in response to delayed IFN by simulating with and without a fixed delay for IFN production from infected cells. Our results suggest that delaying type I IFN production by 5 days yields a roughly 10-fold increase in tissue damage with tissue remaining on day 3 decreasing from 3.9 × 107 cells/ml to 6.5 × 106 cells/ml (Fig 3B), caused by the increase in infected cells and subsequent lack of resistant cells. IFN dynamics were within the observed ranges of systemic IFN- α concentrations from clinical cohorts of hospitalized COVID-19 patients [51] (S7A Fig). As patient IFN-α measurements are only taken after hospitalization and in plasma, there are currently no known sources of early SARS-CoV-2 infection IFN dynamics in humans.
A) Submodel (Eqs 10–16) with all non-IFN cytokines and immune cell interactions set to zero and only considering interactions between virus (V), type I IFN, and susceptible (S), infected (I), resistant (R), and dead (D) epithelial cells. B) Predictions from the simplified model without delayed IFN production (solid lines) versus with a constant delay (τF = 5 days) (dotted lines). Grey circles (left panel): viral loads from SARS-CoV-2 infection in humans in Singapore [48] and Germany [49], digitized from Goyal et al. [50] overlayed with predicted viral dynamics. C-D) Predicted dynamics of infected and dead cells, and unbound and bound IFN concentrations from the simplified model without delayed IFN production (solid lines) versus with a constant delay (τF = 5 days) (dotted lines).
Immunologic determinants of mild and severe disease
Next, to investigate the mechanisms that differentiate mild versus severe disease, we simulated the full model (Eqs. S1-S22) using two different parameter sets. Mild disease dynamics were recreated using estimated parameter values (S1 Table) such that the virus decay rate (dV) and the infected cell death rate (dI) were recalculated to ensure that the maximum death rate of the virus and infected cells did not exceed the value obtained from the reduced viral dynamics model fit (Fig 2). Simulating mild disease, we predicted that all cell populations and cytokines rapidly return to homeostasis, with the immune response effectively clearing virus within 10 days (Fig 4 and S8 Fig).
Mild disease (solid lines) dynamics obtained by using baseline parameter estimates (S1 Table) while severe disease dynamics (dashed lines) were obtained by decreasing the production rate of type I IFN (PF,I) and increasing the production of monocytes (pM,I) and their differentiation to macrophages (). A) Viral load and lung cells concentrations (susceptible, resistant, infected, and dead cells). Solid black line with error bars indicates human data [50] (see Fig 2). B) Immune cell concentrations (inflammatory macrophages, monocytes, neutrophils, and CD8+ T cells). C) Unbound cytokine concentrations (free IL-6, GM-CSF, G-CSF, and type I IFN). Time evolution of all model variables is shown in S8 Fig (including bound cytokine and alveolar macrophages).
Because severe SARS-Cov-2 infection results in lower levels of IFN [51] and increased monocytes [52], we recapitulated severe disease by modulating model parameters relating to these processes, i.e., the rates of IFN production from infected cells, pF,I, and macrophages, , were decreased, and the rate of monocyte recruitment from the bone marrow by infected cells, pM,I, was increased. With these changes, the model predicted a dramatic shift in disease response characterized by a cytokine storm (elevated IL-6, GM-CSF and G-CSF), high ratios of innate to adaptive immune cells, and a marked reduction in healthy viable lung tissue (Fig 4A), whereas changes in viral load remained relatively consistent with mild disease.
In addition, there was a notable increase in the number of inflammatory macrophages (Fig 4B), IL-6, GM-CSF and, importantly, a delayed and reduced IFN peak (Fig 4C). In comparison to mild disease, inflammatory macrophages and neutrophils remained elevated for at least 30 days after initial infection (Fig 4B). Comparing mild and severe disease highlighted differences in the area under the curve (AUC) of macrophages (6 × 104 cells/ml versus 3 × 1011 cells/ml) and neutrophils (2 × 108 cells/ml versus 3 × 1013 cells/ml) over 30 days. Interestingly, inflammation remained high in the severe disease scenario despite the virus being cleared slightly faster (~1 day) than in the case of mild disease (Fig 4A). Further, the peak of inflammatory macrophages increased from ~104 cells/ml to ~106 cells/ml in severe scenarios compared to mild scenarios. The model also accurately predicted that CD8+ T cell dynamics were lower in severe cases, which is indicative of lymphopenia and similar to clinical observations from patients with severe COVID-19 [14,23]. Despite varying only three parameters (pF,I, pM,I, and ηF,MΦ) to generate disparate dynamics, the immune cell and cytokine dynamics were qualitatively in line with clinical observations for IFN-α [51], IL-6 [51,53], and G-CSF [30] (S7B–S7F Fig).
Macrophages, CD8+ T cells, IFN and IL-6 regulate response to SARS-CoV-2 infection
To further understand how the host immune system regulates the response to SARS-CoV-2 infection, we conducted a local sensitivity analysis by varying each parameter individually by ±20% and comparing a set of metrics (see Materials and methods) chosen to provide a comprehensive understanding of each parameter’s impact on the host-pathogen dynamics (S9 Fig). This analysis identified 17 sensitive parameters (Fig 5A) relating to virus productivity (p, δV,N, β, ϵF,I), CD8+ T cell induced epithelial cell apoptosis (δI,T), macrophages, monocyte and CD8+ T cell production (, pM,I, pT,I), IL-6 (
), G-CSF (
), and IFN (pF,I, pF,MΦ,
). Viral load, uninfected cells (tissue), unbound IL-6 and unbound IFN dynamics after 20% decreases in parameters values were found to be similar to those in the original mild disease simulation (Fig 5B–5E).
A local sensitivity analysis was performed by varying each parameter ±20% from its originally estimated value (mild disease parameters in Fig 4) and simulating the model. Predictions were then compared to baseline considering: Maximum viral load (max(V)), maximum concentration of dead cells (max(D)), minimum uninfected live cells (min(S+R)), maximum concentration of inflammatory macrophages (max(MΦI)), maximum number of CD8+ T cells (max(T)), maximum concentration of IL-6 (max(LU)), maximum concentration of type I IFN (max(FU)) and the total exposure to type I IFN (FU exposure). A) Heat map shows the magnitude of the change of each metric from a 20% decrease in the parameter value compared to baseline (i.e. model simulation with no change in the parameter values), where blue signifies the maximum value observed in the output metric and red signifies the minimum value observed (i.e. maximum decrease in that metric). The most sensitive parameters are shown here, for the complete parameter sensitivity results, see S9 Fig. The explicit value of the maximum increase and maximum decrease of each metric is given in the table below. B)-E) time-series dynamics of viral load, tissue (uninfected cells), and unbound IL-6 and IFN given 20% decreases in the noted parameters. Colours of the curves correspond to the colouring of the heatmap in A. Maximal(minimal) concentrations, as in A, are noted in grey boxes. E) is coloured according to IFN exposure. Base: original mild parameters (Fig 4).
Decreasing the rate of IL-6-induced monocyte differentiation into inflammatory macrophages (pMΦ,L) increased the peak concentration of both IL-6 and IFN (Fig 5D and 5E). Notably, changes to parameters that increased the bound IFN concentration, i.e. increasing the binding and production rates ( and pF,I) and decreasing the internalization and clearance rates (
and
) induced significant changes in most metrics (Fig 5A), in particular decreases in
caused an increase in the minimum tissue and IFN (Fig 5C and 5E). The duration of extensive tissue damage (>80% damaged) increased with IFN potency (ϵF,I) (S9 Fig). Less significant changes in the maximum dead cells and minimum tissue can also be seen in the viral infectivity (β) and IL-6 binding kinetic (
and
) parameters. See S9 Fig for complete sensitivity analysis results.
Changes to viral infectivity can be considered linked to two separate, but indistinguishable effects in our model, namely different densities of cellular ACE2 receptors between individuals and changes to viral infectivity through mutations. This would translate to a higher virus infectivity (β) in the model. To predict how increases to viral infectivity (β) values would impact disease trajectories, we varied (β) between 0% and 50% from its original estimated value and assessed predicted outcomes. Although we observed some heterogeneity in response, overall we found that increased infectivity alone is not sufficient to induce significant changes to disease severity (Fig 6), as tissue, cellular, and cytokine concentrations were not dramatically altered from the base case.
A range of viral infectivity rates (β) from 0% (base) to 50% increase were simulated. All other parameters were fixed to their value in S1 Table. The resulting model dynamics for viral load, inflammatory macrophages, CD8+ T cells, unbound IL-6 and unbound IFN were compared, and no significant changes in kinetics were predicted.
In silico knockout investigations predict the impact of immune status on disease outcomes
To better understand the individual contributions of innate immune cell subsets to disease trajectories, we performed a set of in silico experiments in which neutrophil, monocyte, or macrophages were completely removed through knockout (Fig 7). For each scenario, we set the initial cell population to zero (i.e. either N0 = 0, MΦ0 = 0, or M0 = 0) and also fixed the relevant production rates to be zero (i.e. for the neutrophil knockout ; for the macrophage knockout
; and for the monocyte knockout
). In the case of neutrophils, the model predicts that completely removing neutrophils leads to large increases in viral loads, IL-6 concentrations, while depressing macrophages. Similarly, monocyte knockout was predicted to result in significant increases to peak IL-6 concentrations compared to mild disease, whereas macrophage knockout resulted in the collapse of neutrophil numbers and extremely suppressed IL-6 concentrations. Overall, the previously found characteristics of disease severity, including low numbers of unaffected lung tissue cells and large populations of inflammatory macrophages, were not exhibited in any of the knockdown in silico simulations. Given that there is already systemic dysregulation in the severe case, we chose to focus on the mild disease for the knockout simulations to investigate the impact of isolated dysregulation (knockdown and knockout) on immunological trajectories. A comparison of the knockout experiment to severe disease simulations can be found in S10 Fig.
We performed in silico knockout experiments in the mild disease scenario (Fig 4, blue solid lines) by considering complete monocyte knockout (i.e. no monocyte recruitment and M(0) = 0; dark pink dash-dot line), complete macrophage knockout (i.e. not inflammatory macrophage creation via antigen stimulation or monocyte differentiation; light pink dotted line) and complete neutrophil knockout (i.e. no neutrophil recruitment and N(0) = 0; pink dashed line). Dynamics of the in silico knockout are plotted for the A) viral load, B) uninfected cells, C) inflammatory macrophages, D) neutrophils, E) CD8+ T cells relative to uninfected cells and F) unbound IL-6.
Virtual patient cohort identifies heterogeneity in immune dynamics and severity
To better understand the clinical variability in SARS-CoV-2 infection severity [1], we next generated a cohort of 200 virtual patients (see Materials and methods and Fig 8). Here, we make the distinction between a “virtual twin” (i.e., a digital representation of a single clinical patient) and a virtual patient/cohort (multiple patients) generated as explained below. To create each in silico patient, seven patient-specific parameters were sampled from normal distributions with means corresponding to their respective fixed values and standard deviations inferred from clinical observations (Table 1). In doing this, we assumed intrinsic interindividual heterogeneity in monocyte to macrophage differentiation, production of IL-6 by macrophages, recruitment of monocytes by the presence of infected cells, and production of IFN by infected cells, macrophages and monocytes, respectively.
Parameters in the model were first obtained through fitting to data (S1 Table). 1) Parameters relating to macrophage, IL-6 and IFN production (, pL,MΦ, pF,I, pM,I, ηF,MΦ, ϵF,I, and pF,M) were generated from normal distributions with mean equal to their original fitted values and standard deviation informed by experiment observations (see Section S6.1). 2) The model evaluated is then simulated on this parameter set to obtain y(t, p). 3) A simulated annealing algorithm is then used to determine a parameter set that optimises the objective function J(p) (Eq 17). 4) Optimizing the objective function provides a parameter set for which the patient response to SARS-CoV-2 will be within the physiological ranges. This patient is then assigned to the cohort and this process is continued until 200 patients have been generated. Physiological ranges are noted in the bottom box for viral load [37], IFN [51], IL-6 [53] and G-CSF [30].
Seven parameters in the model were deemed patient-specific and were drawn from a normal distribution with mean the parameter value obtained either through fitting or from the literature (S1 Table). The standard deviation (Std Dev) for each normal distribution was informed by values in the literature (see Materials and methods and Supplementary Information Sections S6.1). Initial parameter sampling and new parameters generated through the simulated annealing optimization, were bounded within the interval range noted. All other parameters in the model were fixed to their original value (S1 Table).
Parameters were chosen based on their impact on maximum IL-6 and IFN levels as well as tissue damage observed in the sensitivity analysis (, pL,MΦ, pF,I, pM,I, and ϵF,I; Fig 5). In addition, we designated patient-specific parameters accounting for alternate pathways through which IFN is affected by innate immune cells (ηF,MΦ and pF,M). For the production of IL-6 by macrophages and monocyte to macrophage differentiation via IL-6 stimulation, standard deviations were inferred from IL-6 levels in non-mechanically ventilated patients (mild) and from mechanically ventilated patients (severe) [53] (S7D Fig). Standard deviations for the production of IFN by infected cells were determined from the 95% confidence interval for IFN-α from Trouillet-Assant et al. [51] (S7A and S7B Fig), and, lastly, the standard deviation for the production of IFN by macrophages was obtained from the 95% confidence interval in Sheahan et al. [54] (S7C Fig). The variation in virtual patient responses was then constrained by experimental and clinical viral loads, IFN, neutrophil, IL-6, and G-CSF (Fig 8). The resulting cohort dynamics were within ranges for IFN and IL-6 measurements in asymptomatic to severe COVID-19 patients in the literature [11,17] (S11 Fig).
To quantify disease severity, we introduced an inflammation variable, Ψ, that measured maximum IL-6, neutrophils and tissue damage (Eq 18) and then compared it to individual characteristics of each virtual patient’s disease. We evaluated each virtual patient’s maximum IL-6, CD8+ T cells, and neutrophils; minimum percentage of healthy lung tissue; the time to peak IFN; and total IFN exposure (area under the curve or AUC) within 21 days of infection. Ordering patients by their value of Ψ and plotting the corresponding values for different characteristics showed a clear separation between those with mild disease (Ψ ≤ 3) and those with severe disease (Ψ > 3) (Fig 9A).
200 virtual patients were generated by sampling parameters related to macrophage, IL-6, and IFN production (, pL,MΦ, pF,I, pM,I, ηF,MΦ, ϵF,I, and pF,M) from normal distributions with mean equal to their original values and standard deviation inferred from clinical observations (Fig 8). Each virtual patient had a distinct parameter set optimized to that patient’s dynamics in response to SARS-CoV-2 infection which corresponded to physiological intervals reported in the literature (see Materials and methods). A) Infection and immune response metrics (blue) in individual patients were compared to inflammatory variable Ψ (green). Each point represents an individual patient, ordered according to Ψ. The correlation coefficient (R) and p-value are indicated for each, with α<0.05 denoting significant correlations. B) The effect of exposure dose V0 on maximum IL-6 (a), maximum neutrophil counts (b) and inflammation marker Ψ (c) for V0 = 0.1, 1, 4.5 and 8 log10(copies/ml). In a and b, rows are coloured according to each virtual patient’s inflammation marker value; virtual patients were ordered by the value of Ψ from the baseline scenario in A (V0 = 4.5 log10(copies/ml)). C) Correlations between maximal IFN, IL-6, and T cell concentrations for each patient (circles). Circle colours correspond to the maximal T cell concentration of each patient. D) Parameters most correlated to the IFN peak time were the rates of macrophage production via a) IL-6 (
) and the b) IFN production by infected cells (pF,I). Individual patient values for these parameters are plotted as circles coloured by the patient’s corresponding day of IFN peak (see color bar). Patients were ordered by their inflammation marker Ψ.
In these virtual patients, we investigated predicted disease severity over a range of initial virus exposure doses (V0 = 0.1, 1, 4.5, 8). By comparing outcomes with respect to several immunological biomarkers (e.g. peak IL-6 concentrations, time to IFN peak etc.), we found that virtual patients who were predicted to experience mild COVID-19 in our base case (exposure dose of 4.5 log10 viral copies/mL) experienced mild disease irrespective of the exposure dose size (Fig 9B). In particular, those predisposed to experience severe disease were predicted to experience poor outcomes regardless of inoculation (Fig 9Bc). We observed that patients that had severe disease responses for an inoculation of V0 = 4.5 log10(copies/mL) (i.e. patients 150 to 200) had little to no change in their IFN, IL-6 or neutrophil dynamics under changing inoculation size (Fig 9Ba and 9Bb). However, patients who were prone to have less severe disease responses (i.e. patients 1 to 100) had larger variations in their immune dynamics, in particular IL-6 (Fig 9Ba).
Our model further predicted that patients with higher inflammation had higher IL-6, neutrophil, and inflammatory macrophage concentrations (Fig 9Aa–9Ac), which is somewhat to be expected given that IL-6 is a component of the inflammatory marker Ψ. While the IFN exposure was not significantly stratified by Ψ (Fig 9Ad), the peak of IFN (Fig 9Ae) and CD8+ T cell levels (Fig 9Af) were strongly negatively correlated with the inflammation marker (R = −0.85, p < 1 × 10−9, see Materials and methods). IL-6 was most noticeably correlated with Ψ (R = 0.91, p<1 × 10−9), with a distinct upper bound in the concentration (~80 pg/ml) achieved in 50% of the virtual cohort (Fig 9Aa). There appeared to be a transition phase in inflammation driven by inflammatory macrophage levels where patients with mild inflammation (Ψ < 2.5) had low counts (less than 4 × 104 cells/ml) compared to patients with more severe inflammation (Ψ ≥ 2.5) who had higher levels (p = 1.46 × 10−6; Fig 9Ac). Despite this, patients with moderate inflammation exhibited increased disease markers including delayed IFN peaks and lower CD8+ T cells, compared to patients with mild inflammation (Ψ ≤ 2.5).
A distinct jump in the timing of the IFN peak in the virtual cohort (p <1 × 10−5) was found to be correlated with inflammation (Fig 9Ae), as patients with lower (mild) inflammation (Ψ ≤ 2.5) had peaks at day 2 compared to day 6 in patients with higher (severe) inflammation (Ψ>2.5). Grouping virtual individuals by their time to IFN peak suggests that those with IFN peaks before day 3 of infection also had fewer macrophages (p<1 × 10−5) and larger numbers of CD8+ T cells (p <1 × 10−5). Overall, delays in IFN peak did not cause significant changes to viral load but were sufficient to cause major tissue damage (100x reduction in viable tissue remaining) and over-heightened immune responses (4x increase in maximum IL-6 and GM-CSF concentrations).
Further, examining the relationship between each virtual patient’s maximum IL-6, IFN, and CD8+ T cell concentrations (Fig 9C) identified a weaker correlation between the maximum concentration of CD8+ T cells and IFN (R = 0.24, p = 0.0008) as opposed to with IL-6 (R = −0.86, p < 1 × 10−9). As expected, we also found a positive correlation (R = 0.67, p = 1.58 × 10−8) between the time to peak IFN concentration for each patient and the IFN production rate from infected cells (Fig 9Db). Interestingly, the time to peak IFN for each patient was also strongly related to their rate of IL-6-stimulated monocyte differentiation into macrophages (Fig 9Da). Low IFN production rates were predicted to be the major factor responsible for significantly delayed IFN peaks over 6 days after infection, whereas IFN peaks within 3 days of infection were largely caused by lower rates of monocyte to macrophage differentiation (Fig 9Da).
Discussion
Serial immunological measurements from COVID-19 patients are only beginning to be collected, and the ability to assess initial infection kinetics and the drivers of the ensuing disease course remains limited. The data-driven mechanistic mathematical model and virtual patient cohort developed here is an important platform contributing to investigating immunological drivers of COVID-19. In particular, to recreate severe dynamics, it was sufficient to vary only two processes in the model: the rates of type I IFN production from infected cells and macrophages, and the rate of monocyte recruitment by infected cells. This suggests that the distinction between severe and mild disease may be driven by a limited set of causal regulators that warrant subsequent study. The effect on IFN production may be further exacerbated by autoimmunity against type I IFNs, which has been shown to correlate to life-threatening COVID-19 pneumonia in 2.6% of women and 12.5% of men [18].
Our results show that delaying type I IFN production is sufficient to cause major tissue damage and heightened immune responses, yet it has little impact on peak viral loads. In the severe disease simulation, the viral load was cleared marginally faster (~1 day) in comparison to the mild disease simulation. This finding is supported by recent clinical evidence suggesting that the rate of viral decline may be predictive of disease severity [6,64]. This therefore suggests that viral loads alone may not be a necessary attribute to obtain severe tissue damage. Instead, our model predicts that increases in tissue damage occur through heightened innate immune responses. Further, increases to the viral infectivity rate β, which replicate changes in ACE2 expression between individuals and/or viral mutations causing increases in infectivity, were not predicted to significantly change disease trajectories. This observation is consistent with early SARS-CoV-2 variants, particularly D614G, where a spike protein mutation increased infectivity but did not significantly alter disease outcomes [65] (similar conclusions are emerging for the alpha (B.1.1.7) variant [66,67]). As noted in Table 2, these results suggest future experiments using mutated viral variants, and a separate mechanistic modelling study examining a gradient of ACE2 expression along the respiratory tract and/or within organ systems to understand how within-host variation in ACE2 receptor expression impacts heterogeneity in the same virus. Future work will explore the potential immunopathological effects of newer variants of concern.
Neutrophil and monocyte knockout simulations suggested that excessive suppression of neutrophils early in infection may lead to future hyperinflammation even in otherwise mild COVID-19. The peak ratio between CD8+ T cells and infected cells remained virtually unchanged across the neutrophil, monocyte, and macrophage knockouts we considered, which is indicative of the link between IFN, infected cells, and CD8+ T cells in our model.
Evaluating SARS-CoV-2 infection in a cohort of 200 virtual patients revealed several immunological responses potentially leading to differential disease presentation. Notably, a distinct, emergent switch in the type I IFN response corresponded with late IFN peaks and more severe disease (i.e., higher inflammation Ψ). This supports previous findings that connect a delay in type I IFN with more severe presentations of highly pathogenic coronaviruses infections including SARS-CoV, MERS-CoV, and SARS-CoV-2 [10,13,14], and provides a rational explanation for the finding from a retrospective cohort study that early IFN therapy is associated with better responses [28]. Of note, varying the initial viral inoculum did not impact disease outcomes in virtual patients predicted to have either mild or severe disease. These results seem to suggest that it is an individual’s intrinsic immunological response that dictates the severity of COVID-19. Further, the inflammation marker Ψ can also be applied to patient data, similar to other proposed metrics (e.g. CytoScore [68]).
In our cohort, virtual patients with mild disease tended to achieve peak IFN concentrations approximately 2 days after infection compared to those with severe disease who exhibited higher inflammation and later IFN responses peaking after 5 days. This switch in IFN timing was caused by a 3-fold increase in the rate of monocyte-to-inflammatory macrophage differentiation and decreased production rate of IFN by infected cells. The initial delay of IFN production was caused by increased monocyte-to-macrophage differentiation and this delay was exacerbated by reduced IFN production from infected cells, suggesting that the timing of the IFN peak in a patient may allow for improved stratification into treatment arms designed to target one or both of these responses. The finding that IFN binding was predictive of the duration of lung tissue damage suggests that virus-intrinsic properties and their ability to inhibit receptor mediated binding and endocytosis could delay IFN production and cause downstream increases in IL-6 and GM-CSF, resulting in severe disease. Our results further suggest that lymphopenia is tightly correlated with maximum IL-6 concentration and less dependent on the timing of IFN.
Models by design require simplifying complex dynamics to highlight critical underlying structures, which should then be experimentally or clinically verified. For instance, some unmodeled cytokines have overlapping function and cellular sources to those explicitly modelled here; thus, model predictions may be representative of broader synergistic effects [68] and other cytokines could also be targets of interest. Also, as this work is focused on acute, primary infections, our model does not account for antibody production and we chose not to explore questions related to vaccination efforts nor the potential for cross-reactivity given exposure to other seasonal coronaviruses. Interestingly, recent evidence suggests that antibodies to other seasonal coronaviruses do increase in SARS-CoV-2 but that this does not provide protection [69]. Extending this model in the future to incorporate antibody production and addressing persistent or secondary infection could help shed light on questions about vaccine escape. There are also certain limitations to our virtual cohort simulations. As the true variation of these parameters in humans is unknown and in some cases impossible to measure, we used variation in related data sets for their kinetics. As such, the virtual cohort represents a prognostic guide, and future experiments would be needed to validate the findings in human patients before formal conclusions can be drawn.
Importantly, our approach suggests future avenues of experimental studies through a hypothesis-generation and prediction paradigm, suggesting key mechanisms as promising avenues of investigation (Table 2).
Indeed, the ability of our model to recapitulate severe disease by, in part, regulating monocyte differentiation raises the possibility that patients with low monocyte levels [7] may benefit from treatments that better regulate monocyte differentiation. This is in line with recent studies identifying distinct transcriptional factors as regulators of differentiated monocyte fates in inflammatory conditions [74,75], and clinical observations that monocyte dysregulation is present in severe COVID-19 [76,77]. It also raises the possibility that modulation by exogenous cytokines, including macrophage colony-stimulating factor in combination with IL-4 and tumour necrosis factor-alpha (TNF-α), may be able to direct monocyte differentiation in favour of monocyte-derived dendritic cells and reduce this response [74]. Recently, the neutralization of both TNF-α and IFN-γ has been found to benefit patients with COVID-19 or other cytokine storm-drive syndromes by limiting inflammation and tissue damage [78]. Given that TNF-α also has a secondary benefit on monocyte differentiation, our results support the viability of this avenue of treatment. Caution should be noted, however, given that some previous attempts to regulate host responses by IL-6 blockade have proven unsuccessful [79]. Together, our findings support the idea that early interventions aimed at reducing inflammation are more likely to be beneficial for patients at risk of progressing to severe COVID-19 than attempts to inhibit cytokine storm later in the disease course, given that early IFN responses were found to provoke better controlled immune responses and outcomes in our virtual cohort. It will be essential to characterize both the timing and mechanisms of proposed therapeutic interventions to develop effective treatments to mitigate severe disease.
Materials and methods
Mathematical model of the immune response to SARS-CoV-2
Our model (Eqs. S1-S22) was developed to examine heterogeneous SARS-CoV-2 infection dynamics and explore the immunological drivers of disease severity. The complete model is provided in the Supplementary Information Section S1 along with all variable and parameter descriptions provided in S1 Table, Section S2. Throughout, cytokine and immune cell interactions and effects were described by Hill functions as
(1)
where B is the interacting compound, γ its half-effect value, and h the Hill coefficient [80,81]. Further, for a given cytokine X and cell population Y, the production (recruitment/differentiation) rate of X by Y was denoted by pX,Y and the rate of production of Y by X by pY,X. The half-effect concentration (i.e. γ in Eq 1) of cytokine X on cell population Y was represented by ϵX,Y and the half-effect concentration of cell Y affecting cytokine X was given by ηX,Y. The natural death rate of cell Y was denoted by dY, and the rate of induced death of cell Y by cell Z by δY,Z. Lastly, the carrying capacity concentration of cell Y was denoted by Ymax, and regeneration or proliferation rates by λY.
We modelled virus (V) as produced by infected cells at rate p and cleared via exponential clearance at rate dV, which accounts for all contributions to viral degradation except macrophage- and neutrophil-mediated clearance. Immune-mediated viral clearance via phagocytosis by inflammatory macrophages [82] and neutrophil extracellular traps (NETs—extracellular chromatin fibres produced by neutrophils to control infections) [45,46] was considered to occur at rates and δV,N, respectively. Susceptible epithelial cells (S) grow logistically with per capita proliferation rate λS and carrying capacity Smax, and become infected (I) at rate β. The damage inflicted on epithelial cells by neutrophils was modelled using a Hill function (Eq 1) [81], where neutrophils kill/damage epithelial cells at rate δN through the release of NETs and other antimicrobials proteins [45,46]. The constant ρ (0 < ρ < 1) was included to modulate bystander damage of uninfected cells (S and R).
For the purposes of our study, we only considered type I IFN dynamics (primarily IFN-α, β). Type I IFN (FU and FB) reduces the infectivity and replication capability of viruses by stimulating cells to become resistant to infection [22]. These resistant cells (R) proliferate at a rate equivalent to susceptible cells (λS). The concentration of bound IFN (FB) modulates the creation of infected and resistant cells [19,21,83,84], where increasing the concentration of IFN causes more cells to become resistant to infection and less to become productively infected (I). The potency of this effect is controlled by the half-effect parameter ϵF,I. Following the eclipse phase (which lasts τI hours), productively infected cells (I) were modelled to produce virus before undergoing virus-mediated lysis at rate dI. Although various immune cell subsets contribute to infected cell clearance, we limited our investigation to macrophages and effector CD8+ T cells which induce apoptosis at rates δI,MΦ and δI,T, respectively.
The accumulation of dead cells (D) was assumed to occur through infected cell lysis dI, neutrophil damage/killing of epithelial cells δN, macrophage phagocytosis of infected cells δI,MΦ, macrophage exhaustion δMΦ,D, and CD8+ T cell killing of infected cells δI,T. These dead cells disintegrate relatively quickly [85] at rate dD, and are cleared through phagocytosis by macrophages [86] at rate δD,MΦ.
Resident alveolar macrophages (δΦR) are considered to be replenished at a logistic rate inversely proportion to viral load with maximal rate of λMΦ and half-effect ϵV,MΦ (i.e. as the virus is cleared, the inflammatory macrophage pool replenishes the alveolar macrophage population in the lung). We modelled the transition of alveolar macrophages to inflammatory macrophages (MΦI) as dependent on infected and dead cells, with a maximal rate of aI,MΦ. Resident macrophages die naturally at a rate or due to the clearing of dead cells (exhaustion) [86] at rate δMΦ,D.
Inflammatory macrophages were modelled as produced by three distinct pathways (acting individually or in concert): 1) stimulated tissue-resident macrophages aI,MΦ, (2) GM-CSF-dependent monocyte differentiation, with maximal production pM,G and half effect ϵG,M, and (3) IL-6-dependent monocyte differentiation, with maximal production rate and half-effect
. We assumed that inflammatory macrophages die naturally at rate
or from clearing dead cells at a rate δMΦ,D.
We have previously shown that endogenous cytokine concentrations are not at quasi-equilibrium at homeostasis [87]. Therefore, to describe the pharmacokinetics and pharmacodynamics of cytokine binding and unbinding, we leveraged the framework established in Craig et al. [87] (Fig 1C) for IFN (FB and FU), IL-6 (LB and LU), GM-CSF (GB and GU), and G-CSF (CB and CU). In its general form, this pharmacokinetic relationship is expressed as
(2)
(3)
where YU and YB are free and bound cytokines, Yprod is the rate of endogenous cytokine production, kB and kU are the respective binding and unbinding rates, kint is the internalization rate of bound cytokine, and klin is the elimination rate. Here, POW is a stoichiometric constant, A is a scaling factor and X is the sum of all cells modulated by the cytokine with
(4)
where
is a constant relating the stoichiometry between cytokine molecules and their receptors, K is the number of receptors specific to each cytokine on a cell’s surface and 10n is a factor correcting for cellular units (see Eqs. S19-S22). The molecular weight was calculated in the standard way by dividing the cytokine’s molar mass (MM) by Avogadro’s number (YMW = MM/6.02214 × 1023).
We considered unbound IL-6 (LU) to be produced from productively infected cells, inflammatory macrophages, and monocytes, with bound IL-6 (LB) resulting from binding to receptors on the surface of neutrophils, CD8+ T cells and monocytes. Unbound GM-CSF (GU) was assumed to be produced from inflammatory macrophages and monocytes and bind to receptors on monocytes to create bound GM-CSF (GB). GM-CSF can be produced by CD8+ T cells [88], but this was excluded because it was insignificant to the full system’s dynamics. Unbound G-CSF (CU) is secreted by monocytes, with bound G-CSF (CB) produced via binding to neutrophil receptors. Lastly, because unbound type I IFNs (FU) are known to be produced by multiple cell types in response to viral infection, including lymphocytes, macrophages, endothelial cells and fibroblasts [83], we modelled its unbound production from infected cells, infiltrating/inflammatory macrophages, and monocytes, and its binding to receptors on both CD8+ T cells and infected cells (Fig 1B).
The pharmacokinetics and pharmacodynamics of G-CSF on neutrophils (N) were taken directly from Craig et al. [87]:
(5)
Neutrophil recruitment of bone marrow reservoir neutrophils (NR) was modelled to occur via the bound fraction of G-CSF [89] (CBF = CB(t)/(ACN(t))) at rate which increases towards its maximal value
as a function of increasing G-CSF. During the acute phase of inflammation, endothelial cells produce IL-6 leading to the attraction of neutrophils [90]. This was modelled as recruitment with maximal rate pN,L and half-effect parameter ϵD,L. Neutrophils die at rate dN.
Monocytes (M) are recruited by bound GM-CSF [91], similar to neutrophils (Eq 5), with bone marrow monocytes (MR) recruited at a homeostatic rate . In the presence of GM-CSF, this rate increases towards
. Monocytes are also recruited by the presence of infected cells at a maximal rate of pM,I with half-effect ϵI,M, and subsequently disappear through differentiation into inflammatory macrophages (as above) or death at rate dM.
CD8+ T cells are recruited through antigen presentation on infected cells as a function of infected cell numbers at rate pT,I The constant delay (τT) accounts for the time taken for dendritic cells to activate, migrate to the lymph nodes, activate CD8+ T cells, and the arrival of effector CD8+ T cells at the infection site. CD8+ T cell expansion occurs in response to bound IFN at a maximal rate pT,F with half-effect ϵF,T, and CD8+ T-cell exhaustion occurs with high concentrations of IL-6 [16,17], with half-effect ϵL,T, and apoptosis occurs at rate dT.
Estimating early infection dynamics (‘viral model’)
In an attempt to reduce the degrees of freedom during parameter estimation, we deployed a step-wise approach by isolating subsections of the model. This approach helps mitigate potential issues with parameter identifiability, given that ours is a large, nonlinear model [92–94] and allows us to estimate parameters from multiple data sources [87,95,96]. Other methodologies, including Bayesian computation, are alternative approaches in this context. To begin estimating parameter values from data, we set all immune populations and cytokine concentrations in the full model (Supplementary Information Eqs. S1-S22) to zero (MΦR = MΦI = M = N = T = LU = LB = GU = GB = CU = CB = FU = FB = 0). This gives
(6)
(7)
(8)
(9)
We also assumed there were no resistant cells (R = 0) due to the absence of an IFN equation. This resulted in a simplified ‘viral model’ that considers only virus (V) infection of susceptible cells (S) which creates infected cells (I) after τI days, which the die through lysis, creating dead cells (D).
Type I interferon dynamics during early infection (‘IFN model’)
To study infection dynamics driven uniquely by IFN, we extended Eqs 6–9 by introducing the IFN mechanisms from Eqs. S1-S22, i.e. setting other cytokine and immune cell populations to zero (MΦR = MΦI = M = N = T = LU = LB = GU = GB = CU = CB = 0), giving
(10)
(11)
(12)
(13)
(14)
(15)
(16)
where cells become resistant (R) through IFN (FU and FB). The parameter
was introduced to account for the production of IFN by macrophages and monocytes not explicitly modelled in this reduced system but included in the full system (i.e. pF,M and pF,MΦ in Eq. S17). Previously-fit parameters were then fixed to their estimated values (S1 Table) and the value of
was determined by solving dFU/dt = 0 at homeostasis (i.e. V = I = 0), giving
.
Model calibration and parameter estimation
Model parameters (S1 Table) were obtained either directly from the literature, using the half-life formula (Eq. S23), through fitting effect curves (Eqs. S24-S25) or sub-models (Eqs. S26-S56) to in vitro, in vivo, and clinical data, or by calculating the value that ensured that homeostasis was maintained (Eqs. S57-S70) in the absence of infection. All fitting procedures were performed using MATLAB 2019b functions fmincon or lsqnonlin [97]. Full details are given in the Supplementary Information, and a brief summary is provided below.
Initial concentrations of all unbound cytokines (LU,0, GU,0, CU,0 and FU,0), susceptible cells, resident macrophages, monocytes, neutrophils, and CD8+ T cells (S0, MΦR,0, M0, N0 and T0) were estimated from plasma and lung tissue concentrations in humans [87,98–104] (Section S3.1-S3.2). Parameters for cytokine binding and unbinding kinetics (Eqs 2–4), such as the molecular weight (MM), binding sites per cell (K), binding/unbinding rates (kB and kU), internalization rates for GM-CSF, G-CSF and IFN (kint), and cytokine clearance rates (klin), were estimated both from known values in the literature [84,87,105–115] and previous modelling work [87,116,117] (Section S3.3-S3.5). The stoichiometric constants POW and were both equal to 1 for all cytokines, except for G-CSF for which POW = 1.4608 and
as previously estimated by Craig et al. [87]. Neutrophil and monocyte reservoir dynamics, monocyte differentiation, macrophage activation, and CD8+ T cell recruitment and expansion parameters were primarily estimated from previous mathematical modelling studies [87,118] as well as known values [55,58,81,119] (Section S3.6-S3.8). Immune cell death rates were taken directly from the literature [59,95,120] or estimated from recorded half-lives [60,85,86,121,122] using Eq. S23 (Section S3.9).
To estimate the rates of virus production, decay, infectivity, and infected cell lysis (p, dV, β and dI respectively) in early infections, we initially fit Eqs 6–9 to viral load measurements from SARS-CoV-2 infection in macaques [47] where eight adult rhesus macaques inoculated with 4 × 105 TCID50/ml (3 × 108 genome copies/ml) SARS-CoV-2 [47] (S6 Fig). These parameter values then informed our estimations from hospitalized individuals (Fig 2). We used two data sets of SARS-CoV-2 shedding in the absence of effective treatment from patients in Singapore [40] (n = 3) and patients in Germany [49] (n = 5). In Singapore, samples were obtained with nasopharyngeal swabs whereas viral loads were measured directly from sputum in Germany.
Given the heterogeneity in viral loads from human patients (since available human SARS-CoV-2 viral load data is generally measured from the day of onset of symptoms or after hospitalization without knowing the viral exposure size), Goyal et al. [37] estimated the lag between initial inoculation and first viral measurement for each patient. We used their estimates of the lag and viral production rate p, as well as the estimates for dV, β, and dI from fitting the macaque date (S6 Fig) to estimate the viral parameters dV, β and dI for the human SARS-CoV-2 viral load data. Viral loads below 2 log10(copy/ml) were assumed to be negligible [37]. Estimated parameters for viral decay (dV) and cell lysis (dI) were used as an upper bound for parameter values in the full model to account for additional viral clearance and cell killing of the immune system.
A subset of parameters was obtained through fitting sigmoidal effect curves (Eqs. S24-S25) curves to in vitro and in vivo experiments. These include the IFN inhibition of viral infection and replication [54] (ϵF,I; Section S4.1.1), the half-effect neutrophil concentration for epithelial cell damage [123] (IC50,N; Section S4.1.2), and the half-effect concentrations for monocyte production and differentiation through GM-CSF signalling [124] (ϵG,M and ; Section S4.1.3) see S1 Fig. Other parameters obtained through effect curves were the half-effects for IL-6 production by monocytes [125] and the effect of IL-6 on monocyte differentiation [44] (ηL,M and ϵL,M; Section S4.1.4), and the half-effect of IFN on CD8+ T cell [61] (ϵF,T; Section S4.1.5) and IL-6 on CD8+ T cell expansion [126] (ϵL,T; Section S4.1.6; S2 Fig).
These parameters were then fixed, and remaining parameters were estimated by fitting time-dependent sub-models of Eqs. S1-S22 to relevant data. The proliferation rate of epithelial cells (λS; Section S4.2.1), the internalization rate of IL-6 (; Section S4.2.2), and the rate of neutrophil induced damage (δN; Section S4.2.3) were fit to corresponding time-series measurements [127–129] using exponential rate terms (S2 Fig). Clearance and phagocytosis of infected cells and extracellular virus by inflammatory macrophages (δI,MΦ and δV,MΦ; Section S4.2.4-S4.2.5) were fit to in vitro experiments [86,130] (S2 Fig). Production of IFN by macrophages (pF,MΦ; Section S4.2.6) was obtained by fitting to data measuring IFN-α production [62] (S3 Fig). The parameters regulating the rate of the resident macrophage pool replenishment (λMΦ and ϵV,MΦ; Section S4.2.7) were estimated from our in vivo observations of resident macrophages during influenza virus infection (S3 Fig). GM-CSF production by monocytes (pG,M; S3 Fig, Section S4.2.8), IFN production by infected cells (pF,I; Section S4.2.9), and IL-6 production by infected cells and macrophages (pL,I and pL,MΦ; Section S4.2.10-S4.2.11) were all obtained from fitting reduced versions of Eqs. S1-S22 to in vitro experiments [56,57,131,132] (S4 Fig).
Lastly, any remaining parameters values were obtained by ensuring that homeostasis was maintained in absence of infection (S5 Fig; Section S5). Parameters calculated from homeostasis include the half-effect monocyte concentration for G-CSF production (ηC,M), the production rate of IL-6 and GM-CSF by inflammatory macrophages (pL,MΦ and pG,MΦ), the production rate of monocytes by GM-CSF (pM,G), and the half-effect inflammatory macrophage concentration for IFN production (ηF,MΦ). For some parameters it was not possible to obtain an estimation from the literature, and for these we either set their value equal to an already estimated parameter (), or qualitatively estimated it (ϵI,M, ρ, see S1 Table).
For the ‘IFN model’ (Eqs 10–16), parameters related to virus (p, dV, β and dI), epithelial cell proliferation (λS and Smax), and IFN ( and ϵF,I) were fixed to those in S1 Table.
Numerical simulations
All ODE models were solved using ode45 in MATLAB, and delay differentiation equations (i.e. Eqs. S1-S22) were solved using ddesd in MATLAB.
Sensitivity analysis
We performed a local sensitivity analysis for the full model (Eqs. S1-S22) by individually varying each parameter by ±20% from its estimated value and quantifying the effect on the model’s output. This change was recorded and used to evaluate different metrics representing the inflammatory response to SARS-CoV-2, namely maximum viral load, maximum number of dead cells, minimum uninfected tissue, maximum number of inflammatory macrophages, maximum number of CD8+ T cells, maximum unbound IL-6, maximum unbound IFN, the total exposure (AUC) to type I IFN, number of days the percent of damaged tissue was >80%, and time of unbound type I IFN peak. We quantified the fraction of undamaged tissue by (S + R)/Smax. For each parameter simulation, we recorded the value of the different metrics listed (i.e. maximum viral load etc.). For each metric, we then determined the maximum increase and decrease for that metric, and assigned the grid point for that parameter value a colour based on a linear grid of possible values for that metric.
Virtual patient generation
We next created a virtual cohort of patients to extend the sensitivity analysis in Fig 5 and further interrogate on the causes driving responses for the most sensitive parameters (particularly certain IFN, IL-6, and immune cell related parameters). To generate a cohort of 200 virtual patients, we followed frequently-used quantitative systems pharmacology techniques similar to those of Allen et al. [32] and our previous studies [94,96] wherein individual virtual patients were created by sampling a parameter set p from parameter distributions then simulating the model to verify that each individual’s trajectory was realistic. A subset of parameters (, pL,MΦ, pF,I, pM,I, ηF,MΦ, ϵF,I, and pF,M) was designated as patient-specific after considering the results of the sensitivity analysis (Fig 5 and S9 Fig) and standard deviations inferred from clinical observations (Supplementary Information Section 6.1). Patient-specific parameters were selected using the mvnrand function in Matlab, which samples from a multivariate normal distribution. To avoid the inclusion of unrealistic dynamics, patient parameter sets were then optimized using simulated annealing using the simulannealbnd function in Matlab to ensure predictions fell within physiological ranges for viral load [47], IL-6 [6,53], IFN-α [51], and G-CSF [30] (Fig 8). Posterior distributions from this generation procedure are provided in S10 Fig.
The upper ui and lower li bounds for V, LU, FU and CU were based off these physiological ranges from Munster et al. [47] (viral loads), Herold et al. [53] (IL-6 concentrations), Trouillet-Assant et al. [51] (IFN dynamics), and Liu et al. [7] (G-CSF concentrations) as described in Supplementary Information Section S.6.1. Intervals for each patient-specific parameter set were restricted to four standard deviations from the mean or zero if the lower bound was negative. Given an initial patient specific parameter set p, we used simulated annealing to minimize J(p), i.e.
(17)
where Mi(p) is the model output i evaluated at parameter set p corresponding to the upper and lower bound li and ui (Fig 8).
To quantify disease severity for each patient, we introduced an inflammation variable (Ψ) to account for the combined changes in IL-6 (LU), neutrophils (N), and damaged tissue (S + R), each normalized by the virtual cohort’s average. In this way, Ψ measures an individual’s relative change from the cohort’s baseline, and quantifies the contributions of IL-6, neutrophils, and tissue damage on comparable scales. For a given patient j, the inflammation marker is given by
(18)
where n is the total number of patients in the cohort, and
, and Rj are the unbound IL-6, neutrophils, and susceptible and resistant epithelial cell count, respectively. The threshold value for severe diseases (Ψ = 3) was determined given the distinct “jump” in the delay in peak concentrations (Fig 9).
Statistical analyses
The Pearson correlation coefficient (R) was used to measure the degree of interaction between two variables, with a significance level of α < 0.05 indicating rejection of the hypothesis that there is no relationship between the observed variables. In addition, we used two-sample two-sided t-tests (number of patients < 40) and z-tests (number of patients ≥ 40) at the α < 0.05 significance level to test the hypothesis that there were no differences between sample means.
Supporting information
S1 Fig. Effects of neutrophils on lung epithelial cells, GM-CSF on monocyte production and differentiation, the relationships between monocytes and CD4+ T cells with IL-6, and the influence of IFN on T cell expansion.
A) Using the measurements by Knaapen et al. [123], the inhibitory effect curve E (Eq. S25) was fit to the cell viability of RLE cells under various concentrations of H2O2. B) The stimulatory effect curve E (Eq. S24) was fit to the dose response measurements of blood monoculture cells (3 × 103 cells/dish) with various concentrations of murine recombinant GM-CSF (IU/ml) [124]. C) The stimulatory effect curve E (Eq. S24) was fit to measurements for the monocytic myeloid cell count as a function of GM-CSF [133] D) Eq. S27 fit to time course data of IL-6 production from monocytes [125]. E) IL-6 stimulation of monocyte differentiation to macrophages modelled by the inhibitory effect curve E (Eq. S24) fit to the percentage of CD14+ cells (macrophages) as a function of the number of fibroblasts measured by Chomarat et al. [44]. F) Stimulatory effect curve E (Eq. S24) for IFN-γ stimulation on CD8+ T cells fit to measurements of the signalling in CD8+T cells for varying doses of IFN-γ [61]. Data (black) is plotted as either circles (D & E) or mean and standard deviation error bars (A-C&F); solid blue line: corresponding fit.
https://doi.org/10.1371/journal.ppat.1009753.s002
(TIF)
S2 Fig. Dynamics of IL-6 on T cell expansion, epithelial cell growth, IL-6 internalization, neutrophil-induced damage, and macrophage phagocytosis.
A) Effect curve (Eq. S24) for the IL-6 effect on T cell expansion fit to measurements CD4+ T cells from dilutions of IL-6 by Holsti and Raulet [126]. B) Exponential growth curve fit to the growth of A549 cells [127] C) The internalization rate of IL-6 (Eq. S30) fit to the fraction of internalized IL-6 [128]. D) Exponential decay fit to cell viability after H2O2 administration [129]. E) The macrophage clearance of apoptotic material (Eqs. S31-S33) was fit to the percentage of macrophages that had engulfed material over 25 hours [86]. F) The phagocytosis rate of extracellular virus by macrophages was obtained by fitting Eqs. S34-S35 to the uptake of virus by macrophages measured by Rigden et al. [130]. Data (black) is plotted as either circles (A & F) or mean and standard deviation error bars (B-E); solid blue line: corresponding fit.
https://doi.org/10.1371/journal.ppat.1009753.s003
(TIF)
S3 Fig. Monocyte expansion and type I IFN production by monocytes, alveolar macrophage replenishment after viral infection, and GM-CSF production by monocytes.
A) Eq. S37 fit to time course of proliferation of monocytes in culture [62]. B) Fit of Eqs. S38-S39 to the production of IFN-α by monocytes after 24 hours with RSV as a function of the number of days of pre-culturing (1, 2, 4 or 7) [63]. C) Correlation between infectious virus titre and RT-PCR copy number for influenza A and B measured by Laurie et al. [134] The relative TCID50 compared to the RNA copies is plotted for each virus strain and the mean as a black dashed line. D-E) Fit of Eqs. S40-S42 to viral loads [135] and alveolar macrophages from experimental influenza infections. F) The production of GM-CSF from stimulated monocytes was recorded by Lee et al. [131] Using a simplified version of the full model (Eqs. S43-S46), we obtained the production rates for monocytes and GM-CSF. Data (black) is plotted as either circles/stars (B&F) or mean and standard deviation error bars (A, D-E); solid blue line: corresponding fit.
https://doi.org/10.1371/journal.ppat.1009753.s004
(TIF)
S4 Fig. Production of IFN and IL-6 by infected cells and macrophages.
A) Concentration of IFN-β released by alveolar epithelial cells in response to stimulation with influenza virus recorded at 8, 16 and 24 hours [57]. B-C) IL-6 production by infected cells in response to A) H5NA and B) H7N9, measured by Ye et al. [132] Data (black) is plotted as mean and standard deviation error bars with the corresponding fit (Eqs. S51-S54) in solid blue. D) IL-6 production by macrophages (Eq. S56) in response to stimulation with LPS of varying dosage sizes. Shibata et al. [56] measured the production of IL-6 for different dosages of LPS and fitting the production rate to this data to obtain pL,MΦ, ηL,MΦ.
https://doi.org/10.1371/journal.ppat.1009753.s005
(TIF)
S5 Fig. Homeostatic disease-free system regulation.
A) To confirm that parameters in the model represented realistic immunocompetent individuals in the disease-free scenario, Eqs. S1-S22 were simulated where V0 = 0 and parameters were given by the homeostasis Eqs. S57-S70. The initial concentration of G-CSF was perturbed and compared to simulations of the model at homeostasis. Simulations at homeostasis are represented by solid lines (purple) and perturbed simulations as dashed lines (pink). B) The maximum residual between variables and their initial conditions at day 50 was measured to confirm that the system was stable for perturbations in all immune cells and cytokines.
https://doi.org/10.1371/journal.ppat.1009753.s006
(TIF)
S6 Fig. Viral dynamics model fit to macaque viral data from Munster et al.
[47] A reduced version of the full model (all cytokine and immune cells set to 0, Eqs 6–9) was fit to data from macaques [47] to estimate preliminary viral kinetic parameters. A) Virus (V) infects susceptible cells (S) making infected epithelial cells (I) which then die to produce dead cells (D) and new virus. B) Comparison of predicted viral dynamics compared to observations from 6 animals, with susceptible cell kinetics (left) with predictions of infected and dead cells over time (right). We estimated β, p, dI, V0 and dV from the reduced model in A) fit to data from Munster et al. [47] measuring the viral load in macaques after challenge with SARS-CoV-2.
https://doi.org/10.1371/journal.ppat.1009753.s007
(TIF)
S7 Fig. Model validation against human cytokine and immune cell measurements during SARS-CoV-2 infection.
A) IFN dynamics of the reduced model (Fig 3 Main Text) overlaid with patient IFN-α2 plasma concentrations from Trouillet-Assant et al. [100] The solid line (purple) represents the unbound IFN dynamics from the reduced model (Eqs. 27–33). Individual patient IFN-α2 measurements are plotted as grey circles for IFN-positive patients (n = 21), patients with no IFN measurements (IFN-negative; n = 5) have not been plotted. Healthy volunteer IFN-α2 concentrations are indicated by a grey area. B-F) Mild (solid line) and severe (dashed line) dynamics (Eqs. 27–33 corresponding to simulations in Fig 4 Main Text and S8 Fig compared to corresponding measurements in humans. B-C) Plasma IFN-α and IL-6 in COVID-19 critically ill patients (n = 26) obtained by Trouillet-Assant et al. [100] overlaid with mild and severe unbound IFN (FU(t)) and mild and severe unbound IL-6 (LU(t)) IFN-negative patients (yellow stars) had no IFN-α measurements and IFN-positive patients (grey points) had non-zero IFN-α measurements. Healthy volunteer concentrations are indicated by a grey area. D) IL-6 levels in patients requiring (“Yes”) and not requiring mechanical (“No”) ventilation obtained by Herold et al. [53] overlaid with mild and severe unbound IL-6 dynamics. E) IL-6 concentration in Moderate (“M”) and severe (“S”) COVID-19 patients obtained by Lucas et al. [6]. F) G-CSF plasma concentration obtained by Long et al. [30] in symptomatic “S” and asymptomatic “AS” COVID-19 patients overlaid with corresponding mild and severe model dynamics. G-I) Neutrophils, monocytes and CD8+ T cells in moderate and severe COVID-19 patients normalized by health care worker (HCW) baseline measurements obtained by Lucas et al. [6]. Violin plots are given for the measurements plotted in D-I.
https://doi.org/10.1371/journal.ppat.1009753.s008
(TIF)
S8 Fig. Predicting mild and severe COVID-19 dynamics (all model variables).
Extension of results of mild and severe disease dynamics in Fig 4 Main Text. Mild disease (solid lines) dynamics obtained by using baseline parameter estimates (S1 Table) while severe disease dynamics (dashed lines) were obtained by decreasing the production rate of type I IFN, pF,I, and increasing the production of monocytes, pM,I, and their differentiation to macrophages, ηF,MΦ. A) Lung cells concentrations (susceptible cells S(t), resistant cells R(t), infected cells I(t), dead cells D(t) and virus V(t)). Solid black line with error bars indicates macaque data (see Fig 2 Main Text). B) Immune cell concentrations (resident macrophages MΦR(t), inflammatory macrophages MΦI(t), monocytes M(t), neutrophils N(t) and T cells T(t)). C) Bound and unbound cytokine concentrations (IL-6 unbound LU(t) and bound LB(t), GM-CSF unbound GU(t) and bound GB(t), G-CSF unbound CU(t) and bound CB(t), type I IFN unbound FU(t) and bound FB(t)).
https://doi.org/10.1371/journal.ppat.1009753.s009
(TIF)
S9 Fig. Full analysis of parameters driving COVID-19 severity.
A local sensitivity analysis was performed by varying each parameter ±20% from its originally estimated value and simulating the model. Predictions were then compared to baseline considering: Maximum viral load (max(V)), maximum concentration of dead cells (max(D)), minimum uninfected live cells (min(S+R)), maximum concentration of inflammatory macrophages (max(MΦI)), maximum number of CD8+ T cells (max(T)), maximum concentration of IL-6 (max(LU)), maximum concentration of type I IFN (max(FU)), the total exposure to type I IFN (FU exposure), the number of days damaged tissue was >80% (time (S + R)/Smax)<0.2), and the day type I IFN reached its maximum (day max(FU)). The heatmaps show the fold change of each metric, where blue signifies the minimum value observed and red signifies the maximum value observed, or by the number of days, where light to dark pink signifying increasing number of days from zero. The most sensitive parameters are shown in Fig 5 in the Main Text.
https://doi.org/10.1371/journal.ppat.1009753.s010
(TIF)
S10 Fig. Effects of neutrophil, monocyte, and macrophage knockout on mild disease courses compared with severe disease dynamics.
We performed in silico knockout experiments in the mild disease scenario (Fig 4; solid black line) by considering complete monocyte knockout (i.e. no monocyte recruitment and M(0) = 0; dark pink dash-dot line), complete macrophage knockout (i.e. not inflammatory macrophage creation via antigen stimulation or monocyte differentiation; light pink dotted line) and complete neutrophil knockout (i.e. no neutrophil recruitment and N(0) = 0; pink dashed line). Blue solid lines correspond to mild diseases courses; black solid lines: severe disease. Dynamics of the in silico knockout are plotted for the A) viral load, B) uninfected cells, C) inflammatory macrophages, D) neutrophils, E) CD8+ T cells relative to uninfected cells and F) unbound IL-6. This figure is an extension of Fig 7 in the Main Text.
https://doi.org/10.1371/journal.ppat.1009753.s011
(TIF)
S11 Fig. Parameter distributions of virtual cohort parameters.
Virtual patients were generated by sampling from normal distributions for a subset of model parameters and optimizing parameters using simulated annealing to confirm realistic disease trajectories (see Fig 4). Resulting distributions for each parameter (purple histograms) are shown compared to estimated values for an average patient (dotted black vertical line). The average of the virtual cohort is also plotted (dashed purple line).
https://doi.org/10.1371/journal.ppat.1009753.s012
(TIF)
S12 Fig. Cohort dynamics compared to observed clinical physiological ranges.
Virtual patients were generated so that viral load, IFN, G-CSF and IL-6 concentration were within physiological ranges obtained in the literature. The physiological ranges (Fig 9) of A) human SARS-CoV-2 viral loads [37], B) IL-6 concentrations from patients [53] plotted against the virtual cohort dynamics (grey lines).
https://doi.org/10.1371/journal.ppat.1009753.s013
(TIF)
S1 Table. Parameter values used in the Main Text.
Parameters have been grouped into: (a-e) cell related, (f-k) cytokine related parameters, and (l) initial conditions. Relevant references are given for parameters directly estimated from the literature (“Direct estimate”). Parameters obtained through fitting to data in the literature have the appropriate figure noted and data referenced (“Fit Fig, data”). Remaining parameters were estimated from homeostasis calculation (“Homeostasis”) or qualitatively estimated (“Estimated”). Parameters whose value was taken from another parameters estimation has that parameter noted. Viral load is reported as log(virion copies/ml) and cells have been noted in 109 cells/ml. Time t is in days. The final sub-table (m) is a list of the variables in the model.
https://doi.org/10.1371/journal.ppat.1009753.s014
(PDF)
Acknowledgments
All authors would like to thank Paul Macklin and Thomas Hillen for encouraging collaboration.
References
- 1. Mehta P, McAuley DF, Brown M, Sanchez E, Tattersall RS, Manson JJ. COVID-19: consider cytokine storm syndromes and immunosuppression. Lancet. 2020;395: 1033–1034. pmid:32192578
- 2. Zhou P, Lou Yang X, Wang XG, Hu B, Zhang L, Zhang W, et al. A pneumonia outbreak associated with a new coronavirus of probable bat origin. Nature. 2020;579: 270–273. pmid:32015507
- 3. Qin C, Zhou L, Hu Z, Zhang S, Yang S, Tao Y, et al. Dysregulation of immune response in patients with coronavirus 2019 (COVID-19) in Wuhan, China. Clin Infect Dis. 2020;71: 762–768. pmid:32161940
- 4. Wu F, Zhao S, Yu B, Chen YM, Wang W, Song ZG, et al. A new coronavirus associated with human respiratory disease in China. Nature. 2020;579: 265–269. pmid:32015508
- 5. Stebbing J, Phelan A, Griffin I, Tucker C, Oechsle O, Smith D, et al. COVID-19: combining antiviral and anti-inflammatory treatments. Lancet Infect Dis. 2020;20: 400–402. pmid:32113509
- 6. Lucas C, Wong P, Klein J, Castro TBR, Silva J, Sundaram M, et al. Longitudinal analyses reveal immunological misfiring in severe COVID-19. Nature. 2020;584: 463–469. pmid:32717743
- 7. Liu J, Li S, Liu J, Liang B, Wang X, Wang H, et al. Longitudinal characteristics of lymphocyte responses and cytokine profiles in the peripheral blood of SARS-CoV-2 infected patients. EBioMedicine. 2020;55: 102763. pmid:32361250
- 8. Duployez N, Demonchy J, Berthon C, Goutay J, Caplan M, Moreau AS, et al. Clinico-biological features and clonal hematopoiesis in patients with severe covid-19. Cancers (Basel). 2020;12: 1–11. pmid:32708264
- 9. Nathan C. Neutrophils and COVID-19: Nots, NETs, and knots. J Exp Med. 2020;217. pmid:32678431
- 10. Jamilloux Y, Henry T, Belot A, Viel S, Fauter M, El Jammal T, et al. Should we stimulate or suppress immune responses in COVID-19? Cytokine and anti-cytokine interventions. Autoimmun Rev. 2020;19. pmid:32376392
- 11. Mathew D, Giles JR, Baxter AE, Oldridge DA, Greenplate AR, Wu JE, et al. Deep immune profiling of COVID-19 patients reveals distinct immunotypes with therapeutic implications. Science. 2020;369: eabc8511. pmid:32669297
- 12. Blanco-Melo D, Nilsson-Payant BE, Liu W-C, Uhl S, Hoagland D, Møller R, et al. Imbalanced host response to SARS-CoV-2 drives development of COVID-19. Cell. 2020;181: 1036–1045.e9. pmid:32416070
- 13. Channappanavar R, Fehr AR, Zheng J, Wohlford-Lenane C, Abrahante JE, Mack M, et al. IFN-I response timing relative to virus replication determines MERS coronavirus infection outcomes. J Clin Invest. 2019;129: 3625–3639. pmid:31355779
- 14. Channappanavar R, Fehr AR, Vijay R, Mack M, Zhao J, Meyerholz DK, et al. Dysregulated Type I Interferon and Inflammatory Monocyte-Macrophage Responses Cause Lethal Pneumonia in SARS-CoV-Infected Mice. Cell Host Microbe. 2016;19: 181–193. pmid:26867177
- 15. Li B, Jones LL, Geiger TL. IL-6 Promotes T Cell Proliferation and Expansion under Inflammatory Conditions in Association with Low-Level RORγt Expression. J Immunol. 2018;201: 2934–2946. pmid:30315140
- 16. Liu T, Zhang JJ, Yang Y, Ma H, Li Z, Zhang JJ, et al. The role of interleukin-6 in monitoring severe case of coronavirus disease 2019. EMBO Mol Med. 2020;12: e12421. pmid:32428990
- 17. Laing AG, Lorenc A, del Molino del Barrio I, Das A, Fish M, Monin L, et al. A dynamic COVID-19 immune signature includes associations with poor prognosis. Nat Med. 2020;26: 1623–1635. pmid:32807934
- 18. Bastard P, Rosen LB, Zhang Q, Michailidis E, Hoffmann H-H, Zhang Y, et al. Auto-antibodies against type I IFNs in patients with life-threatening COVID-19. Science. 2020;370: eabd4585. pmid:32972996
- 19. Seo Y-J, Hahm B. Type I interferon modulates the battle of host immune system against viruses. Adv Appl Microbiol. 2010;73: 83–101. pmid:20800760
- 20. Arnaud P. The interferons: pharmacology, mechanism of action, tolerance and side effects. La Rev Med Interne. 2002;23: 449s–458s. pmid:12481400
- 21. Smith AM, Perelson AS. Influenza A virus infection kinetics: quantitative data and models. Wiley Interdiscip Rev Syst Biol Med. 2011;3: 429–445. pmid:21197654
- 22. Park A, Iwasaki A. Type I and Type III Interferons—Induction, Signaling, Evasion, and Application to Combat COVID-19. Cell Host Microbe. 2020;27: 870–878. pmid:32464097
- 23. Sa Ribero M, Jouvenet N, Dreux M, Nisole S. Interplay between SARS-CoV-2 and the type I interferon response. Stapleford K, editor. PLOS Pathog. 2020;16: e1008737. pmid:32726355
- 24. Lokugamage KG, Hage A, de Vries M, Valero-Jimenez AM, Schindewolf C, Dittmann M, et al. Type I Interferon Susceptibility Distinguishes SARS-CoV-2 from SARS-CoV. Williams BRG, editor. J Virol. 2020;94. pmid:32938761
- 25. Mantlo E, Bukreyeva N, Maruyama J, Paessler S, Huang C. Antiviral activities of type I interferons to SARS-CoV-2 infection. Antiviral Res. 2020;179: 104811. pmid:32360182
- 26. Lei C, Qian K, Li T, Zhang S, Fu W, Ding M, et al. Neutralization of SARS-CoV-2 spike pseudotyped virus by recombinant ACE2-Ig. Nat Commun. 2020;11: 2070. pmid:32332765
- 27. Davoudi-Monfared E, Rahmani H, Khalili H, Hajiabdolbaghi M, Salehi M, Abbasian L, et al. A Randomized Clinical Trial of the Efficacy and Safety of Interferon β-1a in Treatment of Severe COVID-19. Antimicrob Agents Chemother. 2020;64. pmid:32661006
- 28. Wang N, Zhan Y, Zhu L, Hou Z, Liu F, Song P, et al. Retrospective Multicenter Cohort Study Shows Early Interferon Therapy Is Associated with Favorable Clinical Responses in COVID-19 Patients. Cell Host Microbe. 2020;28: 455–464.e2. pmid:32707096
- 29. Widagdo W, Okba NMA, Stalin Raj V, Haagmans BL. MERS-coronavirus: From discovery to intervention. One Heal. 2017;3: 11–16. pmid:28616497
- 30. Long Q-X, Tang X-J, Shi Q-L, Li Q, Deng H-J, Yuan J, et al. Clinical and immunological assessment of asymptomatic SARS-CoV-2 infections. Nat Med. 2020;26: 1200–1204. pmid:32555424
- 31. Jenner AL, Aogo RA, Davis CL, Smith AM, Craig M. Leveraging Computational Modelling to Understand Infectious Diseases. Curr Pathobiol Rep. 2020;In press.
- 32. Allen RJ, Rieger TR, Musante CJ. Efficient Generation and Selection of Virtual Populations in Quantitative Systems Pharmacology Models. CPT Pharmacometrics Syst Pharmacol. 2016;5: 140–146. pmid:27069777
- 33. Cassidy T, Humphries AR, Craig M, Mackey MC. Characterizing Chemotherapy-Induced Neutropenia and Monocytopenia Through Mathematical Modelling. Bull Math Biol. 2020;82: 104. pmid:32737602
- 34. Alfonso S, Jenner AL, Craig M. Translational approaches to treating dynamical diseases through in silico clinical trials. Chaos. 2020;30: 123128. pmid:33380031
- 35. Smith AM. Host-pathogen kinetics during influenza infection and coinfection: insights from predictive modeling. Immunol Rev. 2018;285: 97–112. pmid:30129197
- 36. Jenner AL, Yun CO, Yoon A, Coster ACF, Kim PS. Modelling combined virotherapy and immunotherapy: strengthening the antitumour immune response mediated by IL-12 and GM-CSF expression. Lett Biomath. 2018;5: S99–S116.
- 37. Goyal A, Cardozo-Ojeda EF, Schiffer JT. Potency and timing of antiviral therapy as determinants of duration of SARS-CoV-2 shedding and intensity of inflammatory response. Sci Adv. 2020;6: eabc7112. pmid:33097472
- 38. Waghmare A, Krantz EM, Baral S, Vasquez E, Loeffelholz T, Chung EL, et al. Reliability of Self-Sampling for Accurate Assessment of Respiratory Virus Viral and Immunologic Kinetics. J Infect Dis. 2020; 2020.04.03.20051706. pmid:32710762
- 39. Kim KS, Ejima K, Ito Y, Iwanami S, Ohashi H, Koizumi Y, et al. Modelling SARS-CoV-2 dynamics: implications for therapy. medRxiv. 2020; 2020.03.23.20040493.
- 40. Ejima K, Kim KS, Ito Y, Iwanami S, Ohashi H, Koizumi Y, et al. Inferring Timing of Infection Using Within-host SARS-CoV-2 Infection Dynamics Model: Are “Imported Cases” Truly Imported? medRxiv. 2020; 2020.03.30.20040519.
- 41. Sahoo S, Hari K, Jhunjhunwala S, Jolly MK. Mechanistic modeling of the SARS-CoV-2 and immune system interplay unravels design principles for diverse clinicopathological outcomes. bioRxiv. 2020; 2020.05.16.097238.
- 42. Sego T, Aponte-Serrano JO, Gianlupi JF, Heaps SR, Breithaupt K, Brusch L, et al. A modular framework for multiscale multicellular spatial modeling of viral infection, immune response and drug therapy timing and efficacy in epithelial tissues: A multiscale model of viral infection in epithelial tissues. BioRxiv. 2020. pmid:32511367
- 43. Sadler AJ, Williams BRG. Interferon-inducible antiviral effectors. Nat Rev Immunol. 2008;8: 559–568. pmid:18575461
- 44. Chomarat P, Banchereau J, Davoust J, Palucka AK. IL-6 switches the differentiation of monocytes from dendritic cells to macrophages. Nat Immunol. 2000;1: 510–514. pmid:11101873
- 45. Drescher B, Bai F. Neutrophil in viral infections, friend or foe? Virus Res. 2013;171: 1–7. pmid:23178588
- 46. Galani IE, Andreakos E. Neutrophils in viral infections: current concepts and caveats. J Leukoc Biol. 2015;98: 557–564. pmid:26160849
- 47. Munster VJ, Feldmann F, Williamson BN, van Doremalen N, Pérez-Pérez L, Schulz J, et al. Respiratory disease in rhesus macaques inoculated with SARS-CoV-2. Nature. 2020;585: 268–272. pmid:32396922
- 48. Young BE, Ong SWX, Kalimuddin S, Low JG, Tan SY, Loh J, et al. Epidemiologic Features and Clinical Course of Patients Infected with SARS-CoV-2 in Singapore. JAMA—J Am Med Assoc. 2020;323: 1488–1494. pmid:32125362
- 49. Wölfel R, Corman VM, Guggemos W, Seilmaier M, Zange S, Müller MA, et al. Virological assessment of hospitalized patients with COVID-2019. Nature. 2020;581: 465–469. pmid:32235945
- 50. Goyal A, Cardozo-Ojeda EF, Schiffer JT. Potency and timing of antiviral therapy as determinants of duration of SARS CoV-2 shedding and intensity of inflammatory response. medRxiv. 2020; 2020.04.10.20061325. pmid:33097472
- 51. Trouillet-Assant S, Viel S, Gaymard A, Pons S, Richard JC, Perret M, et al. Type I IFN immunoprofiling in COVID-19 patients. J Allergy Clin Immunol. 2020;146: 206–208.e2. pmid:32360285
- 52. Zhao Y, Qin L, Zhang P, Li K, Liang L, Sun J, et al. Longitudinal COVID-19 profiling associates IL-1RA and IL-10 with disease severity and RANTES with mild disease. JCI Insight. 2020;5. pmid:32501293
- 53. Herold T, Jurinovic V, Arnreich C, Lipworth BJ, Hellmuth JC, von Bergwelt-Baildon M, et al. Elevated levels of IL-6 and CRP predict the need for mechanical ventilation in COVID-19. J Allergy Clin Immunol. 2020;146: 128–136.e4. pmid:32425269
- 54. Sheahan TP, Sims AC, Leist SR, Schäfer A, Won J, Brown AJ, et al. Comparative therapeutic efficacy of remdesivir and combination lopinavir, ritonavir, and interferon beta against MERS-CoV. Nat Commun. 2020;11: 222. pmid:31924756
- 55. Kratofil RM, Kubes P, Deniset JF. Monocyte conversion during inflammation and injury. Arterioscler Thromb Vasc Biol. 2017;37: 35–42. pmid:27765768
- 56. Shibata Y, Berclaz PY, Chroneos ZC, Yoshida M, Whitsett JA, Trapnell BC. GM-CSF regulates alveolar macrophage differentiation and innate immunity in the lung through PU.1. Immunity. 2001;15: 557–567. pmid:11672538
- 57. Ioannidis I, Ye F, McNally B, Willette M, Flano E. Toll-Like Receptor Expression and Induction of Type I and Type III Interferons in Primary Airway Epithelial Cells. J Virol. 2013;87: 3261–3270. pmid:23302870
- 58. Pawelek KA, Dor D, Salmeron C, Handel A. Within-host models of high and low pathogenic influenza virus infections: The role of macrophages. PLoS One. 2016;11: 2016. pmid:26918620
- 59. Eftimie R, Eftimie G. Tumour-associated macrophages and oncolytic virotherapies: a mathematical investigation into a complex dynamics. Lett Biomath. 2018;5: S6–S35.
- 60. Patel AA, Zhang Y, Fullerton JN, Boelen L, Rongvaux A, Maini AA, et al. The fate and lifespan of human monocyte subsets in steady state and systemic inflammation. J Exp Med. 2017;214: 1913–1923. pmid:28606987
- 61. Krummel MF, Mahale JN, Uhl LFK, Hardison EA, Mujal AM, Mazet JM, et al. Paracrine costimulation of IFN-γ signaling by integrins modulates CD8 T cell differentiation. Proc Natl Acad Sci U S A. 2018;115: 11585–11590. pmid:30348790
- 62. Ohta M, Okabe T, Ozawa K, Urabe A, Takaku F. 1α,25-Dihydroxy vitamin D3 (calcitriol) stimulates proliferation of human circulating monocytes in vitro. FEBS Lett. 1985;185: 9–13. pmid:3838944
- 63. Krilov LR, Hendry RM, Godfrey E, McIntosh K. Respiratory virus infection of peripheral blood monocytes: Correlation with ageing of cells and interferon production in vitro. J Gen Virol. 1987;68: 1749–1753. pmid:3035067
- 64. Kissler SM, Fauver JR, Mack C, Olesen SW, Tai C, Shiue KY, et al. SARS-CoV-2 viral dynamics in acute infections. medRxiv. 2020.
- 65. Zhang L, Jackson CB, Mou H, Ojha A, Peng H, Quinlan BD, et al. SARS-CoV-2 spike-protein D614G mutation increases virion spike density and infectivity. Nat Commun. 2020;11: 6013. pmid:33243994
- 66. Graham MS, Sudre CH, May A, Antonelli M, Murray B, Varsavsky T, et al. Changes in symptomatology, reinfection, and transmissibility associated with the SARS-CoV-2 variant B.1.1.7: an ecological study. Lancet Public Heal. 2021. pmid:33857453
- 67. Frampton D, Rampling T, Cross A, Bailey H, Heaney J, Byott M, et al. Genomic characteristics and clinical effect of the emergent SARS-CoV-2 B.1.1.7 lineage in London, UK: a whole-genome sequencing and hospital-based cohort study. Lancet Infect Dis. 2021. pmid:33857406
- 68. Brunet-Ratnasingham E, Anand SP, Gantner P, Moquin-Beaudry G, Dyachenko A, Brassard N, et al. Integrated immunovirological profiling validates plasma SARS-CoV-2 RNA as an early predictor of COVID-19 mortality. medRxiv. 2021; 2021.03.18.21253907.
- 69. Anderson EM, Goodwin EC, Verma A, Arevalo CP, Bolton MJ, Weirick ME, et al. Seasonal human coronavirus antibodies are boosted upon SARS-CoV-2 infection but not associated with protection. Cell. 2021. pmid:33631096
- 70. Plante JA, Liu Y, Liu J, Xia H, Johnson BA, Lokugamage KG, et al. Spike mutation D614G alters SARS-CoV-2 fitness. Nature. 2020. pmid:33106671
- 71. Konno Y, Kimura I, Uriu K, Fukushi M, Irie T, Koyanagi Y, et al. SARS-CoV-2 ORF3b Is a Potent Interferon Antagonist Whose Activity Is Increased by a Naturally Occurring Elongation Variant. Cell Rep. 2020;32: 108185. pmid:32941788
- 72. Zheng Y, Zhuang M-W, Han L, Zhang J, Nan M-L, Zhan P, et al. Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) membrane (M) protein inhibits type I and III interferon production by targeting RIG-I/MDA-5 signaling. Signal Transduct Target Ther. 2020;5: 299. pmid:33372174
- 73. O’Brien TR, Thomas DL, Jackson SS, Prokunina-Olsson L, Donnelly RP, Hartmann R. Weak Induction of Interferon Expression by Severe Acute Respiratory Syndrome Coronavirus 2 Supports Clinical Trials of Interferon-λ to Treat Early Coronavirus Disease 2019. Clin Infect Dis. 2020;71: 1410–1412. pmid:32301957
- 74. Goudot C, Coillard A, Villani AC, Gueguen P, Cros A, Sarkizova S, et al. Aryl Hydrocarbon Receptor Controls Monocyte Differentiation into Dendritic Cells versus Macrophages. Immunity. 2017;47: 582–596.e6. pmid:28930664
- 75. Desalegn G, Pabst O. Inflammation triggers immediate rather than progressive changes in monocyte differentiation in the small intestine. Nat Commun. 2019;10: 1–14.
- 76. Rébillard R-M, Charabati M, Grasmuck C, Filali-Mouhim A, Tastet O, Brassard N, et al. Identification of SARS-CoV-2-specific immune alterations in acutely ill patients. J Clin Invest. 2021. pmid:33635833
- 77. Cillo AR, Somasundaram A, Shan F, Cardello C, Workman CJ, Kitsios GD, et al. Bifurcated monocyte states are predictive of mortality in severe COVID-19. bioRxiv. 2021; 2021.02.10.430499. pmid:33594364
- 78. Karki R, Raj BS, Tuladhar S, Williams EP, Zalduondo L, et al. COVID-19 cytokines and the hyperactive immune response: Synergism of TNF-α and IFN-γ in triggering inflammation, tissue damage, and death. bioRxiv. 2020; 2020.10.29.361048. Available: https://doi.org/10.1101/2020.10.29.361048
- 79. Cheng C, Zhang F. Correspondence on: “Interleukin-6 blockade with sarilumab in severe COVID-19 pneumonia with systemic hyperinflammation—An open-label cohort study” by Della-Torre et al. Ann Rheum Dis. 2020;79: 1277–1285. pmid:32620597
- 80. Santillán M. On the Use of the Hill Functions in Mathematical Models of Gene Regulatory Networks. Math Model Nat Phenom. 2008;3: 85–97.
- 81. Norman PS. Immunobiology: The immune system in health and disease. Journal of Allergy and Clinical Immunology. New York: Garland Pub.; 1995.
- 82. Louten J. Virus Transmission and Epidemiology. Essential Human Virology. Elsevier; 2016. pp. 71–92.
- 83. Nagarajan UM. Induction and function of IFNβ during viral and bacterial infection. Critical Reviews in Immunology. Begel House Inc.; 2011. pp. 459–474. pmid:22321107
- 84. Arnaud P. Different interferons: Pharmacology, pharmacokinetics, proposed mechanisms, safety and side effects. Rev Med Interne. 2002;23: 449S–458S.
- 85. Elmore S. Apoptosis: A Review of Programmed Cell Death. Toxicol Pathol. 2007;35: 495–516. pmid:17562483
- 86. Klöditz K, Fadeel B. Three cell deaths and a funeral: macrophage clearance of cells undergoing distinct modes of cell death. Cell Death Discov. 2019;5: 1–9. pmid:30774993
- 87. Craig M, Humphries AR, Mackey MC. A Mathematical Model of Granulopoiesis Incorporating the Negative Feedback Dynamics and Kinetics of G-CSF/Neutrophil Binding and Internalization. Bull Math Biol. 2016;78: 2304–2357. pmid:27324993
- 88. Shi Y, Liu CH, Roberts AI, Das J, Xu G, Ren G, et al. Granulocyte-macrophage colony-stimulating factor (GM-CSF) and T-cell responses: what we do and don’t know. Cell Res. 2006;16: 126–133. pmid:16474424
- 89.
Takei Y, Ando H, Tsutsui K. About the editors. Building Youth for the Future: Suicide Prevention Aspects. San Diego: Academic Press; 2019.
- 90. Scheller J, Chalaris A, Schmidt-Arras D, Rose-John S. The pro-and anti-inflammatory properties of the cytokine interleukin-6. Biochim Biophys Acta (BBA)-Molecular Cell Res. 2011;1813: 878–888. pmid:21296109
- 91. Rösler B, Herold S. Lung epithelial GM-CSF improves host defense function and epithelial repair in influenza virus pneumonia—a new therapeutic strategy? Mol Cell Pediatr. 2016;3: 29. pmid:27480877
- 92. Jenner AL, Frascoli F, Yun CO, Kim PS. Optimising hydrogel release profiles for viro-immunotherapy using oncolytic adenovirus expressing IL-12 and GM-CSF with immature dendritic cells. Appl Sci. 2020;10: 2872.
- 93. Kim PS, Crivelli JJ, Choi I-K, Yun C-O, Wares JR. Quantitative impact of immunomodulation versus oncolysis with cytokine-expressing virus therapeutics. Math Biosci Eng. 2015;12: 841–858. pmid:25974336
- 94. Jenner AL, Cassidy T, Belaid K, Bourgeois-Daigneault M-C, Craig M. In silico trials predict that combination strategies for enhancing vesicular stomatitis oncolytic virus are determined by tumour aggressivity. J Immunother Cancer. 2021;9: e001387. pmid:33608375
- 95. Kim PS, Lee PP, Levy D. Modeling regulation mechanisms in the immune system. J Theor Biol. 2007;246: 33–69. pmid:17270220
- 96. Cassidy T, Craig M. Determinants of combination GM-CSF immunotherapy and oncolytic virotherapy success identified through in silico treatment personalization. PLoS Comput Biol. 2019;15. pmid:31774808
- 97.
The Math Works Inc. Matlab 2019a. MathWorks Inc., Natick, MA. Natick, Massachusetts: Mathwork; 2019.
- 98. Lee J, Kim Y, Lim J, Kim M, Han K. G-CSF and GM-CSF concentrations and receptor expression in peripheral blood leukemic cells from patients with chronic myelogenous leukemia. Ann Clin Lab Sci. 2008;38: 331–337. pmid:18988925
- 99. Rodero MP, Decalf J, Bondet V, Hunt D, Rice GI, Werneke S, et al. Detection of interferon alpha protein reveals differential levels and cellular sources in disease. J Exp Med. 2017;214: 1547–1555. pmid:28420733
- 100. Trouillet-Assant S, Viel S, Gaymard A, Pons S, Richard JC, Perret M, et al. Type I IFN immunoprofiling in COVID-19 patients. J Allergy Clin Immunol. 2020; 4–8. pmid:32360285
- 101. Kasperska-Zajac A, Sztylc J, Machura E, Jop G. Plasma IL-6 concentration correlates with clinical disease activity and serum C-reactive protein concentration in chronic urticaria patients. Clin Exp Allergy. 2011;41: 1386–1391. pmid:21645137
- 102. Uppal SS, Verma S, Dhot PS. Normal Values of CD4 and CD8 Lymphocyte Subsets in Healthy Indian Adults and the Effects of Sex, Age, Ethnicity, and Smoking. Cytom Part B—Clin Cytom. 2003;52: 32–36. pmid:12599179
- 103. Crapo JD, Barry BE, Gehr P, Bachofen M, Weibel ER. Cell number and cell characteristics of the normal human lung. Am Rev Respir Dis. 1982;126: 332–337. pmid:7103258
- 104. CHAFFEY N. Alberts, B., Johnson, A., Lewis, J., Raff, M., Roberts, K. and Walter, P. Molecular biology of the cell. 4th edn. Ann Bot. 2003;91: 401–401.
- 105. Tenhumberg S, Waetzig GH, Chalaris A, Rabe B, Seegert D, Scheller J, et al. Structure-guided optimization of the interleukin-6 trans-signaling antagonist sgp130. J Biol Chem. 2008;283: 27200–27207. pmid:18650419
- 106. Mehra R, Storfer-Isser A, Kirchner HL, Johnson N, Jenny N, Tracy RP, et al. Soluble interleukin 6 receptor: A novel marker of moderate to severe sleep-related breathing disorder. Arch Intern Med. 2006;166: 1725–1731. pmid:16983050
- 107. Stagg J, Wu JH, Bouganim N, Galipeau J. Granulocyte-macrophage colony-stimulating factor and interleukin-2 fusion cDNA for cancer gene immunotherapy. Cancer Res. 2004;64: 8795–8799. pmid:15604233
- 108.
Recombinant human interleukin-6 (with HSA). InvivoGen;
- 109.
Razelle K. Granulocyte-macrophage colony-stimulating factor. Pollock RE WRR et al., editor. Holland-Frei Cancer Medicine. Hamilton (ON): BC Decker; 2003.
- 110.
IFN-beta recombinant protein:: Interferon-beta 1a Recombinant Protein. MyBioSource;
- 111. Chiba S, Tojo A, Kitamura T, Urabe A, Miyazono K, Takaku F. Characterization and molecular features of the cell surface receptor for human granulocyte-macrophage colony-stimulating factor. Leukemia. 1990;4: 29–36. pmid:2153263
- 112. Barreda DR, Hanington PC, Belosevic M. Regulation of myeloid development and function by colony stimulating factors. Dev Comp Immunol. 2004;28: 509–554. pmid:15062647
- 113. Branca AA. Interferon receptors. Vitr Cell Dev Biol—Anim. 1988;24: 155–165. pmid:2450859
- 114. Constantinescu SN, Croze E, Wang C, Murti A, Basu L, Mullersman JE, et al. Role of interferon α/β receptor chain 1 in the structure and transmembrane signaling of the interferon α/β receptor complex. Proc Natl Acad Sci U S A. 1994;91: 9602–9606. pmid:7524081
- 115. Snyers L, Fontaine V, Content J. Modulation of Interleukin-6 Receptors in Human Cells. Ann N Y Acad Sci. 1989;557: 388–395. pmid:2786701
- 116. Nicola NA, Peterson L, Hilton DJ, Metcalf D. Cellular processing of murine colony-stimulating factor (Multi-CSF, GM-CSF, G-CSF) receptors by normal hemopoietic cells and cell lines. Growth Factors. 1988;1: 41–49. pmid:2483336
- 117. Mager DE, Jusko WJ. Receptor-mediated pharmacokinetic/pharmacodynamic model of interferon-β 1a in humans. Pharm Res. 2002;19: 1537–1543. pmid:12425473
- 118. Baral S, Antia R, Dixit NM. A dynamical motif comprising the interactions between antigens and CD8 T cells may underlie the outcomes of viral infections. Proc Natl Acad Sci U S A. 2019;116: 17393–17398. pmid:31413198
- 119. Zhang N, Bevan MJ. CD8+ T Cells: Foot Soldiers of the Immune System. Immunity. 2011;35: 161–168. pmid:21867926
- 120. Craig M, Humphries AR, Mackey MC. An upper bound for the half-removal time of neutrophils from circulation. Blood. 2016;128: 1989–1991. pmid:27578816
- 121. Zent CS, Elliott MR. Maxed out macs: physiologic cell clearance as a function of macrophage phagocytic capacity. FEBS J. 2017;284: 1021–1039. pmid:27863012
- 122. Smith PK, Wang SZ, Dowling KD, Forsyth KD. Leucocyte populations in respiratory syncytial virus-induced bronchiolitis. J Paediatr Child Health. 2001;37: 146–151. pmid:11328469
- 123. Knaapen AM, Seiler F, Schilderman PAEL, Nehls P, Bruch J, Schins RPF, et al. Neutrophils cause oxidative DNA damage in alveolar epithelial cells. Free Radic Biol Med. 1999;27: 234–240. pmid:10443941
- 124. Chen BDM, Clark CR, Chou Th. Granulocyte/macrophage colony-stimulating factor stimulates monocyte and tissue macrophage proliferation and enhances their responsiveness to macrophage colony-stimulating factor. Blood. 1988;71: 997–1002. pmid:2833334
- 125. Alderson MR, Tough TW, Ziegler SF, Grabstein KH. Interleukin 7 induces cytokine secretion and tumoricidal activity by human peripheral blood monocytes. J Exp Med. 1991;173: 923–930. pmid:2007858
- 126. Holsti MA, Raulet DH. IL-6 and IL-1 synergize to stimulate IL-2 production and proliferation of peripheral T cells. J Immunol. 1989;143: 2514–2519. pmid:2571638
- 127. Lawal O, Knobel H, Weda H, Bos LD, Nijsen TME, Goodacre R, et al. Volatile organic compound signature from co-culture of lung epithelial cell line with: Pseudomonas aeruginosa. Analyst. 2018;143: 3148–3155. pmid:29878008
- 128. Nesbitt JE, Fuller GM. Dynamics of interleukin-6 internalization and degradation in rat hepatocytes. J Biol Chem. 1992;267: 5739–5742. pmid:1556093
- 129. Yoshioka Y, Kitao T, Kishino T, Yamamuro A, Maeda S. Nitric Oxide Protects Macrophages from Hydrogen Peroxide-Induced Apoptosis by Inducing the Formation of Catalase. J Immunol. 2006;176: 4675–4681. pmid:16585560
- 130. Rigden RC, Carrasco CP, Summerfield A, McCullough KC. Macrophage phagocytosis of foot-and-mouth disease virus may create infectious carriers. Immunology. 2002;106: 537–548. pmid:12153517
- 131. Lee MT, Kaushansky K, Ralph P, Ladner MB. Differential expression of M-CSF, G-CSF, and GM-CSF by human monocytes. J Leukoc Biol. 1990;47: 275–282. pmid:1689760
- 132. Ye S, Lowther S, Stambas J. Inhibition of Reactive Oxygen Species Production Ameliorates Inflammation Induced by Influenza A Viruses via Upregulation of SOCS1 and SOCS3. J Virol. 2015;89: 2672–2683. pmid:25520513
- 133. Sun L, Rautela J, Delconte RB, Souza-Fonseca-Guimaraes F, Carrington EM, Schenk RL, et al. GM-CSF Quantity Has a Selective Effect on Granulocytic vs. Monocytic Myeloid Development and Function. Front Immunol. 2018;9: 1922. pmid:30210491
- 134. Laurie KL, Guarnaccia TA, Carolan LA, Yan AWC, Aban M, Petrie S, et al. Interval between Infections and Viral Hierarchy Are Determinants of Viral Interference Following Influenza Virus Infection in a Ferret Model. J Infect Dis. 2015;212: 1701–1710. pmid:25943206
- 135. Smith AP, Moquin DJ, Bernhauerova V, Smith AM. Influenza virus infection model with density dependence supports biphasic viral decay. Front Microbiol. 2018;9: 1554. pmid:30042759