COVID-19 virtual patient cohort reveals immune mechanisms driving disease outcomes

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 indicate 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 that was mediated by increased IL-6 and reduced type I IFN signalling. Together, these findings identify biomarkers driving the development of severe COVID-19 and support early interventions aimed at reducing inflammation.

Author summary: Understanding of how the immune system responds 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 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 identify 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 approach to studying novel pathogens using intra-host models.

Figure 1.. Immune response to SARS-CoV-2 infection model schematic.

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_ΦI) 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.

Figure 2.. Viral dynamics model fit to macaque viral data from Munster et al.

[41] A reduced version of the full model (all cytokine and immune cells set to 0, Eqs. 6–9) was fit to data from macaques [41] 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, d₁,V₀ and d_v from the reduced model in A) fit to data from Munster et al. [41] measuring the viral load in macaques after challenge with SARS-CoV-2 (Table S1).

Figure 3.. Delayed type I IFN response impacts heavily on tissue survival in reduced model.

A) Submodel (Eqs. 10–16) with all non-IFN cytokines and immune cell interactions set to zero and only considering interactions between virus (V) 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). Solid black (left panel): viral loads from SARS-CoV-2 infection in macaques by Munster et al. [41] is overlayed with predicted viral dynamics.

Figure 4.. Predicting mild and severe COVID-19 dynamics.

Mild disease (solid lines) dynamics obtained by using baseline parameter estimates (Tables S1) while severe disease dynamics (dashed lines) were obtained by decreasing the production rate of type I IFN (p_F,I) and increasing the production of monocytes (p_M,I) and their differentiation to macrophages (η_F,MΦ). A) Viral load and lung cells concentrations (susceptible, resistant, infected, and dead cells). Solid black line with error bars indicates macaque data [41] (see Figure 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 Figure S7 (including bound cytokine and alveolar macrophages).

Figure 5.. 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(L_U)), maximum concentration of type I IFN (max(F_U)), the total exposure to type I IFN (F_U exposure), the number of days damaged tissue was >80% (time (S + R)/S_max), and the day type I IFN reached its maximum (day max(F_U)). The heatmaps show the magnitude 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 here (for complete parameter sensitivity results, see Figure S8).

Figure 6.. Virtual Cohort of SARS-CoV-2 infected patients.

200 virtual patients were generated by sampling parameters related to macrophage, IL-6, and IFN production ($p_{M{\mathrm{\Phi }}_I,L}$, p_L,MΦ, p_F,I,p_M,I, η_F,MΦ, ϵ_F,I, and p_F,M) from normal distributions with mean equal to their original values and standard deviation inferred from clinical observations (Figure 7). Each virtual patient had a distinct parameter set that was optimized to that patient’s dynamics in response to SARS-CoV-2 infection 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) Parameters most correlated to the IFN peak time were the rates of macrophage production via IL-6 ($p_{M{\mathrm{\Phi }}_I,L}$) and the IFN production by infected cells (p_F,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 are ordered by their inflammation marker Ψ. C) Correlations between maximal IFN, IL-6, and T cell concentrations for each patient (circles). Circle colour corresponds to the maximal T cell concentration of each patient.

Figure 7.. Algorithm for generating virtual patients.

Parameters in the model were first obtained through fitting to data (Table S1). 1) Parameters relating to macrophage, IL-6 and IFN production ($p_{M{\mathrm{\Phi }}_I,L}$, p_L,MΦ, p_F,I, p_M,I, η_F,MΦ, ϵ_F,I, and p_F,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 evaluated 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.16). 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 [41], IFN [42], IL-6 [44] and G-CSF [7].