Quantification of Cytokine Storms During Virus Infections

Highly pathogenic virus infections usually trigger cytokine storms, which may have adverse effects on vital organs and result in high mortalities. The two cytokines interleukin (IL)-4 and interferon (IFN)-γ play key roles in the generation and regulation of cytokine storms. However, it is still unclear whether the cytokine with the largest induction amplitude is the same under different virus infections. It is unknown which is the most critical and whether there are any mathematical formulas that can fit the changing rules of cytokines. Three coronaviruses (SARS-CoV, MERS-CoV, and SARS-CoV-2), three influenza viruses (2009H1N1, H5N1 and H7N9), Ebola virus, human immunodeficiency virus, dengue virus, Zika virus, West Nile virus, hepatitis B virus, hepatitis C virus, and enterovirus 71 were included in this analysis. We retrieved the cytokine fold change (FC), viral load, and clearance rate data from these highly pathogenic virus infections in humans and analyzed the correlations among them. Our analysis showed that interferon-inducible protein (IP)-10, IL-6, IL-8 and IL-17 are the most common cytokines with the largest induction amplitudes. Equations were obtained: the maximum induced cytokine (max) FC = IFN-γ FC × (IFN-γ FC/IL-4 FC) (if IFN-γ FC/IL-4 FC > 1); max FC = IL-4 FC (if IFN-γ FC/IL-4 FC < 1). For IFN-γ-inducible infections, 1.30 × log2 (IFN-γ FC) = log10 (viral load) − 2.48 − 2.83 × (clearance rate). The clinical relevance of cytokines and their antagonists is also discussed.


INTRODUCTION
A large number of studies have demonstrated that viral infections can cause multiple complications in patients that result in multi-organ failure. These may be related to hyper-immune responses to the viruses and may have adverse effects on vital organs and lead to high pathogenicity and mortality (1)(2)(3)(4). For example, the H5N1 influenza virus sets off a cytokine storm, including but not limited to interferon (IFN)-b, interleukin (IL)-6, and interferon-inducible protein (IP)-10 (1-4), and SARS patients showed increased amounts of proinflammatory cytokines in serum [e.g., IL-1b, IL-6, IL-12, IFN-g, IP-10, and monocyte chemotactic protein (MCP) -1], which are associated with pulmonary inflammation and extensive lung damage (5). MERS infection was also reported to induce increased levels of proinflammatory cytokines (e.g., IFNg, tumor necrosis factor (TNF)-a, IL-15, and IL-17) (6). However, it is still unclear whether the cytokine with the largest induction amplitude is the same under different virus infections. It is unknown which among these cytokines is the most critical and whether there are any mathematical models that can fit the dynamics of the cytokines.

Data Extraction and Quality Assessment
We conducted a literature search of peer-reviewed publications in the PubMed electronic database from its inception to October 15, 2020. Only human infections by highly pathogenic viruses were included in this analysis. They were three coronaviruses (SARS-CoV, MERS-CoV, and SARS-CoV-2), three influenza viruses (2009H1N1, H5N1 and H7N9), Ebola virus, human immunodeficiency virus, dengue virus, Zika virus, West Nile virus, hepatitis B virus, hepatitis C virus, and enterovirus 71. Low-pathogenic seasonal H3N2 infection was investigated as the control group.
To get cytokine change information, the following three sets of keywords were employed for the literature search: "patient" (keyword 1), "cytokine" (keyword 2) and individual virus name (keyword 3) ( Figure S1). We retrieved the full text of the potentially eligible studies and examined the full-text reports to obtain information about fold changes of cytokines relative to the healthy (or convalescent) control group. The articles without individual cytokine levels in the control group were excluded. Studies where n < 5 (such as some case reports) or three or more FC > 100 (which should be outliers) were also excluded because that the data may be unreliable. If there were multiple sampling time points, the highest (peak) value was selected. For the case where two or more references showed the same cytokines, the data with the largest n value were selected. Some cytokine levels in the control group were extremely low (close to zero) and were therefore excluded from subsequent analysis or adjusted to the normal levels as indicated in the references (Table S1). If none of these articles included either IFN-g or IL-4, then more articles that included these two cytokines were searched (keyword 2 "patient" was replaced with "IL-4" or "IFN-g"). The final set of papers with their cytokine information are listed in Tables S1 and S2.
To get viral load information, the following three sets of keywords were employed for the literature search: "patient" (keyword 1), "viral load" (keyword 2) and individual virus name (keyword 3) ( Figure S2). Different reports presented viral loads in different ways, such as cycle threshold (Ct) values of the RT-PCR analysis, log 10 viral RNA copies, or absolute viral titers. Only the data expressed as RNA copies/mL were recorded. The papers where n < 5 (such as some case reports) were excluded. For respiratory viruses, viral loads in throat swabs, sputum, or lower respiratory tracts were recorded. For the other viruses, serum viral loads were recorded. If there were multiple sampling sites, the highest value was selected; if there were multiple sampling time points, the highest (peak) value was selected. In case two or more references showed the same cytokines, the data with the largest n value were selected (Table S3).
To obtain virus clearance time information, the following three sets of keywords were employed for the literature search: "patient" (keyword 1), "viral shedding or clearance" (keyword 2) and individual virus name (keyword 3) ( Figure S3). The duration of viral shedding or the time from onset of symptoms to negative PCR result (50 th percentiles for the time until the loss of virus RNA detection) for each virus was recorded. The papers where n < 5 (such as some case reports) or only with the time of hospital stay or the duration of fever were excluded. For respiratory viruses, median durations of virus in respiratory samples were recorded. For other viruses, data in serum samples were recorded. In case two or more references showed the viral load or the clearance time from the same infection, the data with the largest n value were selected (Table S3).
We did not consider the deviation effects of drug treatments because for a certain virus infection with a certain symptom, the drug treatment is usually fixed. For example, both oseltamivir and antibiotics are usually given for severe influenza virus infections (independent of hemagglutinin types) at the time of admission (7). Cytokine data, viral load, and clearance rate were usually acquired from the same patients with the same drug treatment.

Data Processing and Statistical Analysis
Correlations between IFN-g FC and max FC, IL-4 FC and max FC, and max FC were calculated.
Correlations between viral load and log 2 (IFN-g FC) or between clearance rate and log 2 (IFN-g FC) were calculated. For Ebola survivors, Ebola fatalities, dengue hemorrhagic fever patients, and HCV patients, the clearance rates were close to zero (Table S3). Under these four conditions, the IFN-g levels may be mainly affected by the viral loads and not by the clearance rates. We calculated the linear regression between IFN-g FC and log 10 (viral load), log 2 (IFN-g FC) and log 10 (viral load), or log 10 (IFN-g FC) and the log 10 (viral load). Log 2 (IFN-g FC) got the best linear regression: log 10 (viral load) -2.48 = 1.30 × log 2 (IFN-g FC).
Given that viral load and clearance rate are positively and negatively correlated with IFN-g FC, respectively, we presume that 1.30 × log 2 (IFN-g FC) = log 10 (viral load) -2.48 − coefficient × clearance rate. Among all the viruses selected in this study, the fastest clearance happens to 2009H1N1 that 6.84 log 10 virus particles (8) could be cleared within 5 days in the mild patients (9) with only a 1.3 fold-induction to IFN-g (10); thus the coefficient could be calculated as 2.83. Then, the correlation between [log 10 (viral load) − 2.48 − 2.83 × (clearance rate)]/1.30 and log 2 (IFN-g FC) was calculated.
The F-test was performed to analyze all the correlations and determine whether the data pairs fit the regression model. The regression equation, the correlation coefficient, and the P-value were obtained by using SPSS v19.0 and Microsoft Excel 2013. Pvalue threshold and R 2 threshold for statistical significance for claims of correlations were 0.05 and 0.5 respectively. Low-pathogenic seasonal H3N2 infection was used as a control sample ( Figures S1-S3 and Tables S1-S4). According to the severity of symptoms, a total of 25 virus infection cases were summarized, and the two cytokines with the most significant increases were recorded in each case (Figures 1, S4-S8). Most of the cytokines increased after viral infections, with a maximum increase of 102 times, and the increase in some individual patients was up to 200 times (e.g., IL-8 in Ebola fatalities; Figure  S6) (11). However, some other cytokines were significantly inhibited by viral infections, such as granulocyte colony-stimulating factor (G-CSF) and vascular endothelial growth factor (VEGF) in Ebola virus infection (Figures 1, S6) (11) and IFN-g and IL-17 in HIV infection (Figures 1, S8) (12). Among the 50 cytokines with the largest induction amplitudes, IP-10, L-6, IL-8, and IL-17 appeared with the highest frequency (5 out of 50). IFN-g and IL-4 appeared four times; and MCP-1 appeared three times (Table S2). It is worth noting that the cytokines with the largest induction amplitudes were not necessarily the same in mild patients and severe patients with the same virus infection. For example, in EV71 patients with encephalitis, IL-13 and IL-4 were the most increased cytokines, whereas in EV71 patients without encephalitis, IL-22 and macrophage inflammatory protein-1a (MIP-1a) were the most increased cytokines (Table S2 and Figure S8) (13,14).

Variations of IFN-g and IL-4 Determine the Maximum Amplitude of All Cytokines
Cytokines produced by T helper (Th) cells are of critical importance for the outcome of many infectious diseases. Producing the "right" set of cytokines in response to an infectious agent can be a matter of life or death. Although the Thl/Th2 dichotomy (mutual antagonistic loop) is an oversimplification, it has proven useful in the analysis of immune responses to infections (Infante-Duart and Kamradt, 1999; Paludan et al., 1998). The two cytokines IL-4 and IFN-g play major roles in the generation and regulation of immune responses. Central in this respect are their mutually antagonistic functions. IFN-g plays a key role in the inhibition of Th2-cell differentiation and Th1-cell stabilization; IL-4 promotes Th2-cell differentiation and stability and inhibits Th1-cell differentiation (15,16). A significant correlation between IFN-g FC and the maximum induced cytokine (max) FC was found (R 2 = 0.697; Figure 2A), and no significant correlation between IL-4 FC and max FC was found (R 2 = 0.228; Figure 2B).

Viral Load and Clearance Rate Are Positively and Negatively Correlated With IFN-g Fold-Changes
Next, the determinants of IFN-g changes were investigated. For the condition that IFN-g increased after the infection, viral load was weakly positive-correlated with log 2 (IFN-g FC) (R 2 = 0.382; Figure 3A), and virus clearance rate was negatively correlated with log 2 (IFN-g FC) (R 2 = 0.678; Figure 3B). Ebola survivors need 158 days to clear the virus, and Ebola fatalities result when the individual cannot clear the virus before death (17). Patients with primary dengue hemorrhagic fever need a long time (far more than 12 days) to clear the viral protein NS1 (18); 85% of HCV patients cannot clear the virus within 9 months (15% of patients cleared HCV spontaneously within 108 days) (19,20). Their clearance rates are close to zero. Thus, under these four conditions, the IFN-g levels may be mainly affected by the viral loads but not by the clearance rates. Then, we calculated the linear regression between IFN-g FC and log 10 (viral load), log 2 (IFN-g FC) and log 10 (viral load), or log 10 (IFN-g FC) and the log 10 (viral load). Log 2 (IFN-g FC) got the best linear regression ( Figure 3C). The following regression equation was obtained: log 10 (viral load) -2.48 = 1.30 × log 2 (IFN-g FC), which means that every 10-fold increase in viral load would result in 1.7-fold induction to IFN-g. The intercept of 2.48 means that when viral load < 2.48 log 10 , IFN-g could not be enhanced by the infection. In fact, viral loads less than 10 2 -10 3 are usually below the PCR detection limit and do not induce cytokine storms (17)(18)(19)(20)(21).
Given that viral load and clearance rate are positively and negatively correlated with IFN-g FC, respectively, we presume that 1.30 × log 2 (IFN-g FC) = log 10 (viral load) -2.48 − coefficient × clearance rate. Among all the viruses selected in this study, the fastest clearance happens to 2009H1N1 that 6.84 log 10 virus particles (8) could be cleared within 5 days in the mild patients (9) with only a 1.3 fold-induction to IFN-g (10); the coefficient could be calculated as 2.83. For the IFN-g-inducible conditions, the correlation between [log 10 (viral load) − 2.48 − 2.83 × (clearance rate)]/1.30 and log 2 (IFN-g FC) was calculated and a high correlation coefficient R 2 = 0.916 was obtained ( Figures 3D, E), which suggests that the inducing amplitude of IFN-g could be predicted according to the viral load and the clearance rate. However, no such correlations could be found for the condition that IFN-g decreased after the infection (such as HIV, HBV, and EV71 infections).

DISCUSSION
The cytokines with the largest induction amplitudes are summarized in this study. Their antagonists may be used clinically to control the inflammatory responses. However, under different virus infections, the most significantly increased cytokines are usually different. There is often a difference between mild patients and severe patients with the same virus infection. Therefore, the kind of cytokine antagonists that should be used to control the inflammatory responses depends on individual situations. Many previous studies focused on IL-6, which may contribute to disease exacerbation, and some therapeutic approaches based on anti-IL-6 biologics have been proposed (22,23) and validated (24). However, our analysis showed that besides IL-6, IP-10, IL-17, and IL-8 are the most common cytokines with the largest induction amplitudes (e.g., IL-8 may be the most significantly induced cytokine in SARS-CoV-2 infection) (25).
Antagonists against IP-10, IL-17, or IL-8 have not attracted enough attention in clinical practice. For the Th1-type infections, IFN-g is the most important determinant of cytokine storm severity. Therefore, IFN-g antagonists may be used as candidate drugs to control the inflammatory response. In most cases of Th1-type infections (IFN-g FC/IL-4 FC > 1), IFN-g is increased. There is one exception: although IFN-g decreased after HBV infection, IL-4 decreased further (26); so, the ratio of IFN-g FC/IL-4 FC was still greater than 1 (Table S2). In this case, IFN-g antagonists should not be used. Given that max FC = (IFN-g FC) 2 /IL-4 FC, enhancement to IL-4 may also prevent the cytokine storm. However, excessive IL-4 may suppress the inflammation too much and generate detrimental effects (15,16). Over-inhibition to IFN-g may also down-regulate immunity seriously, and therefore delay the clearance of the virus (15,16). We should be very cautious when using either IFN-g antagonists or IL-4.
For Th2-type infections (IFN-g FC/IL-4 FC < 1), IFN-g might be an ideal drug. IFN-g has been proved to show some therapeutic effects to HIV patients (27). Our results suggest that IFN-g treatment may also be extended to other Th2-type infections, such as the EV71 infection.
Our analysis suggests that viral replication and clearance rate are the decisive factors in inducing IFN-g and other cytokines, thus determining severity of the cytokine storm ( Figure 4). For example, the viral load of Ebola survivors was two orders of magnitude lower than Ebola fatalities (28), so a difference of 1.7 times (23.3 vs 40) in IFN-g FC and a difference of 2.5 times (41.5 vs 102) in max FC were observed (Table S2) (11). 6.84 log 10 2009H1N1 virus particles could be cleared within 5 days (8,9), whereas clearance of 7.07 log 10 H7N9 virus requires 19.7 days (7,29). Correspondingly, 1.3-fold and 9.3fold inductions to IFN-g were observed, respectively, and a difference of 5 times (2 vs. 10) in max FC was observed (Table  S2) (10,30). The viral load usually reaches a peak at 1-10 days post the infection, and the peak of cytokines often appears at the same time with the peak load or appears in 1-7 days post the peak load (Tables S1-S4). The best way to prevent a cytokine storm is to reduce the viral load or accelerate the virus clearance, which suggests that antiviral therapy should be started as early as possible.
A large number of studies have developed cytokine changes into boolean models, and then, converted the models to ordinary differential equations (ODEs) that helped to reveal the behaviour of the various cytokines (31)(32)(33)(34)(35)(36). However, these reports mainly focused on dynamics of individual cytokines (31,32,(34)(35)(36) or cytokine groups (such as pro-inflammatory cytokine group and anti-inflammatory cytokine group) (33). Interactions among IL-4, IFN-g and cytokines with the largest induction amplitudes have not been modeled before. Nevertheless, our mathematical models are merely based on linear regression and quite simplified. Up till now, there are too few references for some infections (such as H3N2, H5N1, H7N9, Ebola, WNV and EV71) to construct a more sophisticated model. With more extensive research, more quantitative data about cytokine storms would be published, and then more precisely mathematical models may be established.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
SY had the idea for and designed the study and had full access to all data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis. SY, and SCJ contributed to writing of the manuscript. All authors contributed to the literature search, data acquisition, data analysis, statistical analysis and reviewed and approved the final version.

ACKNOWLEDGMENTS
We thank LetPub for its linguistic assistance during the preparation of this manuscript.

SUPPLEMENTARY MATERIAL
The