The C-ImmSim model accounts for differences in the HLA haplotype when determining which peptides are presented by antigen-presenting cells. To this end, for each HLA molecule, it takes in input a list of such peptides together with a propensity of each of them to bind to the HLA. This list is computed by using third-party immunoinformatics tools as described in the next section (Computing the peptide immunogenicity).

The “HLA haplotype freq search” in the “Allele Frequency Net Database” (45) was used to select two HLA-A, two HLA-B and two DRB alleles which are most prevalent in the US population (46). The result pointed to the following alleles: HLA-A*02:01, HLA-A*24:02, HLA-B*35:01, HLA-B*40:02, DRB1*07:01, and DRB1*15:01.

The strain of SARS-CoV-2 used in this study corresponds to the reference sequence NCBI Reference Sequence: NC_045512.2. The primary structure of these proteins has been used to identify cytotoxic T peptides (CTL peptides) and helper T peptides (HTL peptides). To this aim, we have employed two immunoinformatics tools. In particular, for the definition of CTL epitopes, the “ANN 4.0 prediction method” in the online tool MHC-I binding prediction of the IEDB Analysis Resource (47) was used for the prediction of 9-mer long CTL peptides which had an affinity for the chosen set of HLA class I alleles (i.e., HLA-A*02:01, HLA-A*24:02, HLA-B*35:01 and HLA-B*40:02) (48–50). The peptides were classified as strong, moderate, and weak binders based on the peptide percentile rank and IC50 value. Peptides with IC50 values <50 nM were considered to have high affinity, <500 nM intermediate affinity, and <5000 nM low affinity towards a particular HLA allele. Also, the lower the percentile rank, the greater is meant the affinity (48–50). The list of peptides is reported in Appendix B of the Supplementary Materials .

For what concerns the HTL peptides, the NetMHCIIpan 3.2 server (51) was used for the prediction of 9-mer long HTL peptides which had an affinity for the HLA class II alleles (i.e., DRB1*07:01 and DRB1*15:01) used in this study (52). The predicted peptides were classified as strong, intermediate, and non-binders based on the concept of percentile rank as given by NetMHCII pan 3.2 server with a threshold value set at 2, 10, and >10%, respectively. In other words, peptides with percentile rank ≤2 were considered as strong binders whereas a percentile rank between 2 and 10% designate moderate binders; peptides with percentile score >10 are considered to be non-binders (52). The list of CTL and HTL peptides and the relative affinity score is reported in Appendix C of the Supplementary Materials .

It is widely accepted that aging is accompanied by remodeling of the immune system. With time, there is a decline in overall immune efficacy, which manifests itself as an increased vulnerability to infectious diseases, a diminished response to antigens (including vaccines), and a susceptibility to inflammatory diseases. The most important age-associated immune alteration is the reduction in the number of peripheral blood naïve cells, accompanied by a relative increase in the frequency of antigen-specific memory cells. These two alterations are extensively reported in the literature and account for the immune repertoire reduction (53, 54). Along with the process called “inflammaging”, the reduction of immune repertoire is considered the hallmark of immunosenescence (55).

To model the reduction of immune efficacy we first defined the parameter “immunological competence” IC∈(0,1] and assumed it in a simple linear relationship with age. Specifically, we set IC ≡ IC(age) = –α age + 1 with the value of the parameter α = 45 · 10^-4 determined using epidemiological data as described below. Given the age, the parameter IC is then used to modulate both innate and adaptive immunity as follows: i) the phagocytic activity of macrophages and dendritic cells, represented by a probability to capture a viral particle, is rescaled respectively as p_M = IC · u and p_DC = IC · v, where u~U _[a,b] and v~U _{[5× a ,5× b ]} are two random variables uniformly distributed in the ranges [a, b] and [5 · a, 5 · b] with a = 25 × 10^-4 and b = 10^-2; ii) as for the adaptive immunity it is adjusted according to the immunological competence parameter IC by decreasing the lymphocyte counts (hence B, Th, and Tc) to reflect a reduction in the repertoire of “naïve” cells with immunological history due to accumulation of memory cells filling the immunological compartment. In particular, the number of white blood cells N is computed as N ∼ I C · N ( μ , σ 2 ) = N ( I C · μ , I C 2 · σ 2 ) where N ( μ , σ 2 ) is a normal distribution with average μ and standard deviation σ (for each lymphocyte type B, Th and Tc) chosen to reflect the reference leukocyte formula for an average healthy human adult (see Figure 2 ) (56).

The infection and the dynamical features of the SARS-CoV-2 viral strain have been characterized by two parameters: i) V ₀ corresponding to the infectious viral load at time zero, and, ii) the affinity of the virus spike protein to the ACE2 receptor on target cells, called p_A . In particular, V 0 ∼ U [ c 1 , d 1 ] has been taken randomly in the interval [c ₁ = 5, d ₁ = 5 · 10⁵], while p A ∼ U [ c 2 , d 2 ] in the interval [c ₂ = 10^-3, d ₂ = 10^-1].

Upon choice of the age-class determining the immunocompetence value IC(age) hence p_A and p_DC as well as the lymphocyte counts, the simulations depict the immune-virus competition eventually culminating in a successful, or not, virus-clearing response controlling its growth. Sometimes this control is not perfectly efficient. In those cases, the result is a longer viral persistence possibly going much beyond the length of the observation period of 30 days (cf. Figure 4 ).

The sequence of events from viral infection leading to a full-fledged immune response is detailed in Appendix A of the Supplementary Materials . At each time step of the simulation, C-ImmSim dumps all variables allowing for a detailed analysis of the dynamics. A full output example of a simulation is reported and described in Appendix D of the Supplementary Materials .

The analogy of some simulation variable to a realistic clinical endpoint allows stratifying the in-silico patients for a more concrete interpretation of the results. Based on the viral load observed at day 30 (indicated by V ₃₀) and a threshold θ we assume to identify the virtual patients in one of the three following classes:

According to this assumption/definition, by choosing the cutoff θ =120 viral particles per micro-liter of simulated volume, we obtain the stratification of the virtual individuals shown in Figure 3A . Of all in-silico individuals, we get 4.3% of critical cases (broken down in age classes in Figure 3A ), 46.8% of partially recovered ( Figure 3B ), and 48.8% recovered cases ( Figure 3C ). These figures sound very much in line with current epidemiological statistics when considering that the recovered cases here simulated include asymptomatics (57).

This result shows an interesting and surprisingly high fraction of in-silico patients in the partially recovered class. This class includes patients that, at the end of the simulated period of thirty days, are still positive albeit manifesting an active immune response, regardless of being asymptomatic or not. This question is discussed below.

These special cases can be better examined in Figure 4 , which shows four distinct exemplifying runs with different outcomes. In Figure 4A the viremia is shown as a function of time. Red lines correspond to individuals who reach the critical condition V ₃₀ > θ thus falling in the class “critical”. The green line corresponds to a viral clearance corresponding to a fully recovered case, and the blue line shows a situation in which the virus is not completely cleared but stays below the threshold value θ. This case corresponds to one of what we call partially recovered as it represents virtual individuals that produce an immune response (cf. same figure, Figure 4A showing the corresponding antibody titers) which turns out to be insufficient to clear the virus. These “unresolved infections” include asymptomatic cases and are worth the further analysis described below.

To note that the two examples of critical outcomes (red curves) originate from a quite different initial viral load. Also, to note that the fully recovered (green) case starts with a viral load that is higher than one critical case, still the immune response manages to control the infection. The blue curve shows a partially recovered case which greatly decreases the viral load (inset plot of Figure 4A ) but does not clear it completely.

It is worth clarifying that the term “symptom” has no meaning in the in-silico framework unless a link between model variables and possible clinical endpoints is drawn. Also, we should note that we have no concept of comorbidity here that would help in defining the “status” of the virtual patient. To overcome this limitation, besides the viral load at day 30, we think up the following quantities (or variables) as clinical endpoints: (a) the damage in the epithelial compartment, namely, percent of virus-target cells that are dead at the time of observation as a surrogate marker of vascular permeability; (b) the concentration of pyrogenic cytokines as a surrogate marker of fever, i.e., Prostaglandins TNFa, IL-1, and IL-6 causes fever people get varying degree of severity with COVID-19.

Of these two potential surrogate markers of criticality, the first appears more appropriate. While the amount of pyrogenic cytokines (surrogate clinical endpoint b.) correlates with the severity of the disease, the most striking difference in the critical cases versus non-critical (i.e., recovered plus partially recovered) is seen when comparing the accumulated damage in the epithelial compartments (φ) computed as the fraction of depleted epithelial cells due to the immune cytotoxicity of SARS-CoV-2 infected cells during the whole observation period (surrogate clinical endpoint a.),

where E(t) is the epithelial count per microliter of simulated volume and T is the time horizon of 30 days (note that φ ∈ [0,1]). Indeed, when we plot the distribution of φ for the critical and non-critical cases (i.e., partially recovered plus fully recovered) we obtain what is shown in Figure 5 . The plot clearly shows that for critical in-silico patients the damage is much more pronounced than for non-critical ones. This prompt us to use the threshold φ_c = 0.63 to set apart patients who have mild infections [about 80% as in reality (60–62)] to those having a severe disease [15% with dyspnea, hypoxia, lung changes on images (63)] or critical illness [5%, respiratory failure, shock, multi organic dysfunction, cytokine storm syndrome (64)]. Therefore, we label patients with φ ≥ φ_c as symptomatic while those with φ < φ_c asymptomatic. According to this further stratification patients who are still positive (i.e., 0 ≤ V ₃₀ < θ) and have no symptoms (i.e., φ < φ_c ) account for about 44% of the simulations which is in line with current estimates of asymptomatic incidence (65).

We tested the correlation between the antigen abundance (or infective viral load V ₀) and the severity of the infection. The Mann-Whitney-Wilcoxon (MWW) test shows a significant (i.e., p-value< 10^-3) difference between infecting viral load V ₀ in the three classes critical, partially recovered, recovered. In particular, we find that a higher V ₀ is a strong correlate of disease severity (66). Of interest is the fact that there is no significant difference among age groups, that is, V ₀ is not predictive of the disease severity for the age (67) (MWW test p-value>0.05).

Significantly, in most critically ill patients, SARS-CoV-2 infection is associated with a severe clinical inflammatory picture based on a severe cytokine storm that is mainly characterized by elevated plasma concentrations of interleukin 6 (68). In this scenario, it seems that IL-6 owns an important driving role in the cytokine storm, leading to lung damage and reduced survival (69).

The simulation agrees with this finding as reported in Figure 6 . Plotting the peak value of the viral load (i.e., the maximum value attained in the observed period) versus the logarithm of the integral of IL-6 over the whole period (cf. eq(2) in next section 3.5), we see a positive correlation regardless of the outcome.

For all age classes, as shown in Figure 7A , a critical clinical course is associated with a significantly higher concentration of pro-inflammatory cytokine IL-6. Figure 7B shows the same information for all age classes lumped together. The cytokine concentration on the y-axis is calculated as the integral over the whole simulated period (definition in eq(2) of section 3.5). The positive correlation between inflammation and severity of the clinical course is a consequence of the struggle of the immune system to cope with the infection. However what Figure 6 reports is a generic higher production of IL-6 in younger individuals compared to elders. The explanation of this outcome becomes visible following the line of consequences starting from a stronger cytotoxic activity (see Figure 12B below) that killing infected cells cause a stronger release of danger signal to which macrophages respond by secreting IL-6. Since younger have a higher immunological competence (IC), they respond with both stronger cytotoxic response and better innate (i.e., macrophage) activity. The result is the somehow counterintuitive observation that while younger individuals have a higher degree of inflammation, they report a smaller propensity to experience severe outcomes.

What is observed in the production of IL-6 in younger individuals extends to all cytokines produced during the response to the inflammation. We find that, in general, cytokines’ cumulative production during the whole simulated period correlates inversely with the age. Calling c_x the concentration of cytokines x in the simulated volume, where x is one of IL-6, D, IFNg, and IL-12, we define

the cumulative value of cytokines in the whole observation period. Figure 8 shows σ_{IL 6}, σ_D , σ_IFNg , and σ_{IL 12} against age. What Figure 8 shows is that there is a clear reduction of cytokines’ production with age. This, similarly to what was discussed in the previous section, is due to the reduced immune activity which indeed is determined by a reduced immunological competence with the age (71). To justify the apparent contrast with the fact that elder acute infected individuals are more prone to experience a cytokine storm, we should openly regard one of the limitations of the model, namely, the lack of further cytokine feedbacks that are activated during extended pneumonia.

The expression of IFNg by CD4 tends to be lower in severe cases than in moderate cases as shown in Figure 9A . This is in agreement with (72).

The inverse correlation of interferon-gamma (IFNg) with disease severity is observed in all age groups ( Figure 9A and also in Figure 9C when summing all age classes). Interestingly, recovered and partially recovered do not show a meaningful difference when compared to the critical cases ( Figure 9C ).

IFNg is released by natural killer (NK) cells upon bystander stimulation by danger signals (Rule n.5 in Appendix A of Supplementary Materials ) which, in turn, is released by infected/injured epithelial cells upon viral infection (Rule n.3) and when killed by cytotoxic cells (Rule n.18). This means that in young individuals a prompt activation of NK cells due to higher immunological competence, and a stronger cytotoxic response killing infected cells, are sufficient to control the “acute” production of danger signal impacting on the production of IFNg.

If we use the cumulative value of inflammatory cytokines as a variable, namely, σ_IFNγ + σ_{IL 6} + D + σ_TNFa (i.e., the sum of the integrals) and we use the accumulated damage in the epithelial compartments, that is, the fraction of depleted epithelial cells due to the immune cytotoxicity φ defined in eq(1) (cf. section 3.2) as the discriminating criteria between symptomatic and asymptomatic, we observe what is shown in Figure 10 .

Compared with asymptomatic cases, the symptomatic ones more frequently have markedly higher levels of inflammatory cytokines. The difference of the virtual patients in the two classes is statistically significant (MWW test, p-value< 10^-3), which is in line with the clinical finding of higher inflammatory levels in severe disease progressions (72). This result seems to contrast what is stated in section 3.5, namely that younger individuals produce more cytokines but have a less severe course of the disease. However, the explanation provided by the simulation is that those who deal with the infection more rapidly (including asymptomatic) produce, overall, a smaller amount of cytokines and therefore have a lower risk of “cytokine storm”. On the other hand, an inconclusive immune response leads to chronic inflammation, hence pronounced symptoms (e.g., extended epithelial damage) and ultimately in a cytokine storm.

It has been suggested that in younger individuals several factors contribute to the lower numbers of patients observed with severe disease, namely: lower number of ACE receptors, overlapping immunity against coronaviruses, and a more efficient intact immune system.

Figure 11 shows that the immune response is quicker in younger virtual individuals compared to elder ones. The speed of the immune response is calculated in terms of the time [in days) the viral load V(t) reaches its maximum [indicated t V m a x where t:V(t) = V^max and V max = max t V ( t ) ] and starts to decline due to the immune response. Figure 11A shows the distribution of t V m a x for each age class. Younger individuals develop a faster response and consequently, the virus is cleared earlier. This is shown in Figure 11B which plots the distribution of the time (in days) it takes the immune response to decrease the viral load below the threshold θ whenever this happens (the cases for which V(t) > θ,∀t, are not counted in this statistics). Figure 11B is in line with the fact that younger individuals mount a quicker immune response that is generally more efficient than those in elder people thus eradicating the virus in a shorter time.

Figure 12 shows that younger individuals have a higher production of antibodies when compared to elder individuals. This is evident for both critical and recovered (partial or fully recovered) individuals. However, the most striking observation when considering the difference between recovered and critical cases is the gap in antibody titers present in virtually all age classes ( Figure 12A ). This indicates a strong protective role of the humoral response making a split between recovered and critical patients.

Figure 12B shows the corresponding statistics for the cytotoxic T cell (peak value) count per age-class and critical status. This plot consistently evidences that youngers have a stronger response than elders. Interestingly, in contrast to the humoral response, in all age classes, the cytotoxic response in critical individuals is higher than in recovered ones revealing the attempt of the immune system to counterbalance the inefficacy of the humoral response.

Moreover, a further view at the antibody titers reveals that its peak value [i.e., V m a x = max t V ( t ) ] correlates inversely with the clearance time (t_θ ), that is, faster responses are obtained with lower production of antibodies (cf. Figure 13 ). This is in line with the hypothesis that asymptomatic individuals develop a rapid but mild response which clears the infection (73).

It should be noted, however, that there are still substantial uncertainties on the data available due to the variable diagnostic accuracy of different serological tests for COVID-19, therefore more well-designed large clinical studies are warranted to address this matter (74, 75). Interestingly, the same cannot be said when considering peak values of Tc counts, that is, the cytotoxic response does not correlate, either positively or negatively, with time-to-clearance (not shown).

We have used a logistic regression classification to see if by measuring the antibody titers and the CTL counts at day t < 30 we can infer the outcome at the end of the simulated period of t = 30 days. We call V(t = 30) = V ₃₀ the viral load on day 30 after the infection.

Formally, the logistic regression classifier uses the data set {(x_i , y_i )}_{i=1… m} where x i = ( x i ( 1 ) , x i ( 2 ) ) is the feature vector consisting of the normalized cytotoxic T-cell lymphocytes count, x i ( 1 ) = T c ( t ) and the antibody titers at day t, x i ( 2 ) = A b ( t ) . While y_i = 0 if the corresponding run has V ₃₀ ≤ θ and y_i = 1 if the corresponding runs has V ₃₀ > θ.

Figure 14A shows the features x_i corresponding to the recovered cases (i.e., y_i = 0) represented as yellow circles and the critical cases (i.e., y_i = 1) corresponding to black daggers. This panel shows the best separation curve found after training a logistic regression model on the training sample x_i = [Tc(25), Ab(25)], namely, Tc count and antibody titers measured at day 25 from infection.

Figure 14A shows the data set after the classification in recovered and critical and the separation curve. In the figure, the data set corresponds to the observation at day 30 while the analysis has been conducted at different time points. Figure 14B shows how the classification accuracy increases over time. In this panel, we plot the Sørensen-Dice coefficient [most known as the F1 score (76)] which increases when the assessment is made by using features (i.e., Tc and Ab measurements) later as the infection and corresponding immune response develops. Interestingly, before day 10 after the infection it is not possible to find a meaningful classification criterion that predicts the outcome, while the Sørensen-Dice coefficient increases to a high value already at day 25 indicating that 25 days after infection, the level of immune activation represented by the antibody titers and cytotoxic counts is predictive of the clinical outcome.