Understanding the immune response to COVID-19 is challenging due to its high variability among individuals. To identify differentially expressed cytokines between the deteriorating and recovering phases, we analyzed the Electronic Health Records (EHR) and cytokine profile data in a COVID-19 cohort of 444 infected patients and 145 non-infected healthy individuals. We categorized each patient's progression into Deterioration Phase (DP) and Recovery Phase (RP) using longitudinal neutrophil, lymphocyte and lactate dehydrogenase levels. A random forest model was built using healthy and severe patients to compute the contribution of each cytokine toward disease progression using Shapley Additive Explanations (SHAP). SHAP values were used for supervised clustering to identify DP and RP-related samples and their associated cytokines. The identified clusters effectively discriminated DP and RP samples, suggesting that the cytokine profiles differed between deteriorating and recovering health conditions. Especially, CXCL10, GDF15, PTX3, and TNFSF10 were differentially expressed between the DP and RP samples, which are involved in the JAK-STAT, NF- B, and MAPK signaling pathways contributing to the inflammatory response. Collectively, we characterized the immune response in terms of disease progression of COVID-19 with deteriorating and recovering health conditions.