1.

Introduction

Coronavirus disease 2019 (COVID-19) is caused by the severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2). SARS-CoV-2 is absorbed by the mucosae of the respiratory tract and targets bronchial cells and pneumocytes. Once SARS-CoV-2 is internalized into the host cell, it is sensed by receptors especially present in native cells and macrophages, a subtype of immune cell.^1,2 The binding of viral particles to their receptors leads to a pro-inflammatory response by the infected cells. The innate immune response is the immediate mechanism of defense by the immune system and is led by neutrophils and macrophages, which release molecules (cytokines) that act locally to stop the aggressor as well as attract lymphocytes responsible for viral recognition and regulation of the immune cell response.^2,3 This phase, led by both T and B lymphocytes, is known as the adaptive immune response.⁴

In lung tissue from severely ill COVID-19 patients, some cytokines are overexpressed⁴ (cytokine storm), which is thought to result in severe collateral damage in the affected tissues, driving organ failure.² This cell network between immune and surrounding cells has been morphologically explored previously, mainly in routine histopathology postmortem tissue from people who have died on account of COVID-19.^5–13 Additionally, the most common immune cell response consists of a predominance of interstitial T-lymphocytes.¹⁴ In contrast, B cells are consistently few while neutrophils and macrophages are heterogeneously reported.¹⁴

H1N1 is a viral disease that is still causing thousands of deaths.¹⁵ Previous works^15–17 have highlighted the importance of comparing the differences between COVID-19 and H1N1 since they share common etiologies and similar clinical symptoms, occur in the same season, and were considered pandemics in recent times. Given the novelty of COVID-19 and the non-specific nature of clinical presentation of both diseases, there is a need for better understanding how these diseases impact the immune system.

When comparing COVID-19 with H1N1, it has been found that subjects with COVID-19 present (1) higher median age of onset; (2) more prevalence of symptoms such as non-productive cough, fatigue, and gastrointestinal disturbances; (3) more ground-glass opacities in radiology images; (4) worse response to some oxygen therapies; and (5) wider spectrum of pharmacotherapy with different responses.¹⁵ A possible reason for these could be on account of the differentially expressed immune response between H1N1 and SARS-CoV-2.^15,18,19 In the lung, the spectrum of histopathological findings in COVID-19 includes diffuse alveolar damage, fibrin microthrombi, desquamative pneumocytes, and superposed changes of bronchopneumonia; all of them accompanied by T lymphocytic infiltration, in mild to moderate degree. In H1N1, additional findings such as necrosis of main airways, extensive hemorrhage,²⁰ superimposed acute bronchopneumonia, and a considerable amount of interstitial macrophages accompanying lymphocytic infiltration^21,22 have been more consistently found. However, none of these studies used objective metrics for comparing differences in architectural patterns of the immune response in both diseases using pathology slides.

Although visual assessment of hematoxylin–eosin (H&E)-stained tissues and immunohistochemistry of autopsied specimens in both COVID-19 and H1N1^14,23 provide a gestalt interpretation, it involves some degree of subjectivity and does not yield quantitative measurements. Computational pathology, image analysis, and machine learning approaches have been previously employed for quantitative characterization of the immune response in digitized H&E tissue pathology images of diseases such as breast and lung cancer.^24–27

In this preliminary work, we introduce a new computational imaging framework to characterize the disease microenvironment from H&E whole slide images (WSIs) of postmortem lung tissue belonging to nine patients who succumbed to COVID-19. First, image processing and machine learning were employed to automatically identify two nuclei types in each image tile: “lymphocytes” and “non-lymphocytes.” Although non-lymphocytes in lung tissue include pneumocytes, macrophages, neutrophils, among other nucleated cells, in this study, these individual cell types were not differentiated but treated as corresponding to a single family. Then nuclei clusters for both types were built based on the proximity of the individual nuclei. In addition, we present a new set of quantitative metrics that capture the density and spatial interplay between clusters of the same and different nuclei types. Unlike the previous work employing nuclei-graph methods,²⁸ this approach takes into account inter- and intraprimitive relationships (i.e., relationships between nuclei types of the same or of different families), thereby enabling the exploration of complex cellular interactions within the microenvironment. Although, in this work, we showcase the utility of this approach in characterizing the interplay and spatial relationships between two different nuclei types (i.e., lymphocytes and non-lymphocytes), the approach can be extended to other nuclei types and also to a larger number of nuclei families (i.e., $>2$). Additionally, we apply this approach to characterize the architectural nuclei arrangement found in two patients who died due to H1N1 and identify a set of image features relating to the immune morphology that were most different between H1N1 and COVID-19 patients.

2.

Materials and Methods

2.2.

Automatic Identification of Lymphocytes

For automatic identification of lymphocytes (Fig. 1), image color normalization³¹ was applied to compensate staining variations of slides acquired from different institutions [Fig. 1(c)]. Next, the method of Veta et al.,³² which has been used by previous works,^26,33,34 was used for segmenting all image nuclei [Fig. 1(d)]. This method is a watershed-based model specifically tuned for nuclei segmentation of H&E images. In an independent testing set containing 18 histopathology cases, it showed a sensitivity of 0.85 and mean predictive value of 0.89.³²

Fig. 1

Procedure for automatic detection of lymphocytes. Top row: (a) A WSI; (b) an enlarged field of view extracted from the WSI; (c) field of view after applying normalization of color; (d) automatic segmentation of nuclei; and (e) detection of lymphocytes. Bottom row: For each segmented nucleus, a set of features related to texture, color, and shape are extracted. Then a previously trained classifier (a support vector machine) uses such features to label each nucleus as either a lymphocyte or a non-lymphocyte.

Finally, the model by Corredor et al.³⁵ was employed to classify each nucleus as either “lymphocyte” or “non-lymphocyte.” This model consists in a support vector machine with linear kernel that employs image-derived features of each segmented nucleus related to texture, shape, and color (Fig. 1, bottom row). Authors of the method report that in an independent testing set containing 1026 lung nuclei (manually annotated and validated by an expert pathologist), it yielded a precision of 0.89 and recall of 0.83.³⁵

2.3.

Characterizing the Spatial Architecture of Lymphocytes and Surrounding Nucleated Cells

Once lymphocytes and non-lymphocytes were automatically identified, different clusters for each family were built based on the distance between nuclei. First, lymphocytes were connected with each other if the distance between them is below a predefined threshold. A very large threshold value generates just a single complete graph, which limits analyzing interactions within the disease microenvironment; whereas a small value (close to zero) produces multiple sparse graphs, which approximates to analyze individual nuclei. For this reason, in this study, the value of the threshold was empirically determined to be 50 pixels ($\sim 12.5\mu \mathrm{m}$) because it showed reasonable trade-off in the number of generated graphs. This action was carried out for each lymphocyte until all lymphocytes have been visited and potentially connected with corresponding proximal lymphocytes. This results in $n$ disconnected subgraphs of lymphocytes [Fig. 2(b)]. The same process is carried out for non-lymphocytes, thereby generating $m$ disconnected subgraphs of non-lymphocytes. Then a set $P$ of polygons is constructed and defined by $P=\{P_1^L,P_2^L,\dots ,P_n^L,P_1^{NL},P_2^{NL},\dots P_m^{NL}\}$, where $P_i^L$ is the convex hull³⁶ that contains all the nuclei of the $i$’th lymphocyte subgraph and $P_j^{LN}$ is the convex hull of the $j$’th non-lymphocyte subgraph [Fig. 2(c)]; each polygon represents a cluster. The convex hull has been used previously to model the area of influence, for example, the spatial extent of epidemics at the outbreak stage.^37,38 Given the dynamics of the disease microenvironment, in this work, the convex hull concept was employed to model the area of influence of nuclei, taking into account their movements within the tissue, the associated cytoplasm, and connecting stroma.

Fig. 2

Construction of nuclei clusters: (a) identification of two families of nuclei (blue: lymphocytes and green: non-lymphocytes); (b) linking of nuclei based on proximity; and (c) resulting nuclei clusters after drawing the convex hull that contains all the linked nuclei. Bottom cell of (b) illustrates three nuclei $\alpha$, $\beta$, and $\gamma$. Since the distance between $\alpha$ and $\beta$ ($d_1$) is less than the threshold $\theta$, they are linked. On the contrary, the distance between $\alpha$ and $\gamma$ ($d_2$) is larger than $\theta$, so they are not connected. As a result, they belong to different clusters. The nuclei that are not linked to other nuclei cannot conform clusters, so they are not considered for the analysis.

Based on the built clusters, six groups of features were derived.

• 1. Number of clusters of each nuclei type. $N_f$ represents the number of clusters of the $f$’th nuclei type, where $f\in \{1,2,\dots ,F\}$, where $F$ is the total number of families ($F\ge 2$).

• 2. Density of each cluster. $d_k=n_k/A_k$, where $d_k$ is the density of the $k$’th cluster, $n_k$ is the number of nodes constituting the cluster $k$, and $A_k$ is the area of cluster $k$, i.e., the total number of pixels covered by the convex hull $C$ of the cluster nodes.

• 3. Intersection between clusters. $A_i\cap A_j$, $\frac{A_i\cap A_j}{A_i}$, $\frac{A_i\cap A_j}{A_j}$, $\frac{A_i\cap A_j}{(A_i+A_j)/2}$, where $A_i\cap A_j$ is the intersection between clusters $i$ and $j$, $A_i$ is the area of cluster $i$, and $A_j$ is the area of cluster $j$ [see Fig. 3(a)].

• 4. Diversity of clusters surrounding a specific cluster. $S_{f,k}/N_k$, where $S_{f,k}$ represents the numbers of clusters of nuclei type $f$ surrounding the $k\text{-}i$’th cluster ($f\in \{1,\dots ,F\}$) and $S_k$ represents the total number of clusters surrounding cluster $k$. This measure is applied for a predefined number of neighbors [see Fig. 3(b)].

• 5. Measures from global graphs built for each nuclei type. A complete and undirected graph $G=(O,E,W)$ was defined for each nuclei type, where $O=\{o_1,\dots ,o_F\}$ is the set of nodes corresponding to the convex hull centroids of each nuclei type, $E=\{e_1,e_2,\dots ,e_m\}$ is the set of edges connecting such centroids such that $\{(o_i,o_j)\in E:\forall o_i,o_j\in O,i,j\in \{1,\dots ,F\},i\ne j\}$, and $W=\{w_1,w_2,\dots ,w_m\}$ is a set of weights proportional to the length of each $e\in E$. To extract information about the arrangement of cluster centroids of a nuclei type, subgraphs representing the Voronoi diagram $G_V$, Delaunay triangulation $G_D$ [Fig. 3(c)], and MST $G_{\mathrm{MST}}$ were constructed. From such subgraphs, different features are computed, such as area, perimeter, and chord length of Voronoi polygons, area and side length of Delaunay triangles, edge length of MST, among others.

• 6. Compactness of cluster centroids of each nuclei type. Given the graph $G$ previously defined, we computed the compactness of centroids as $c_i=\sum _{j=1}^{N_{f-1}}\frac{1}{D_{i,j}}*b_{\theta }$, where $c_i$ is the compactness measure of the $i$’th node, $N_f$ is the number of nodes of nuclei type $f$, $D_{i,j}$ is the Euclidean distance between nodes $i$ and $j$, $\theta$ is a threshold value, and $b_{\theta }$ is a binary value that is 0 when $D_{i,j}>\theta$ and 1 otherwise [see Fig. 3(d)].

Fig. 3

Computation of some representative features: (a) intersecting area between clusters; (b) cluster neighborhood diversity; (c) Delaunay graph built for a specific nuclei type; and (d) node compactness for a specific nuclei type. The color bar represents the grouping measurement, in which $H$ stands for nodes highly grouped, i.e., very close to multiple nodes, while $L$ stands for sparsely clustered nodes, i.e., isolated or far from other nodes.

For characterizing each patch, a set of statistics (total, mean, standard deviation, median, minimum, maximum, skewness, and kurtosis) were computed for each feature.

3.

Results

3.1.

Experiment 1: Analyzing Spatial Patterns of Immune Response in COVID-19

Figure 4 presents the resulting t-SNE plot for all the patches of COVID-19 patients. No obvious grouping of patient patches was observed, suggesting that spatial architecture features are sparsely distributed and there is not a drastic difference from one patient to another.

Fig. 4

Unsupervised clustering of spatial architectural features of immune clusters of image patches from COVID-19 patients. Each point represents a patch and each color represents a different patient. This plot was generated using the t-SNE algorithm that started out with 350 features and reduced the dimensionality to 2 for facilitating visualization.

Some lymphocyte-related features such as the area and density of lymphocyte clusters, intersection between clusters of lymphocytes non-lymphocytes, graph-based measures, and compactness of clusters were also analyzed using violin plots (Fig. 5). The median area of clusters of lymphocytes is around 3000 pixels ($\sim 0.33\%$ of the patch area), indicating that clusters of lymphocytes are generally small [Fig. 5(a)]. Similarly, the median intersected area between clusters of lymphocytes and non-lymphocytes was $\sim 2500\text{pixels}$ ($\sim 83\%$ of the median area of lymphocyte clusters), suggesting that a considerable area of the lymphocyte clusters was overlapping with non-lymphocyte clusters [Fig. 5(c)]. Finally, graph-based and compactness measures [Figs. 5(d)–5(f)] revealed that lymphocyte clusters tend to be circumscribed by other lymphocyte clusters.

Fig. 5

Violin plots comparing the variation of some lymphocyte-related features across different COVID-19 patients ($N=9$). For the sake of visualization, some extreme outliers were removed from subfigures (a)–(c). The dashed line indicates the mean value over all nine patients. (a) Median area (pixels) of lymphocyte clusters per patch; (b) median density of lymphocyte clusters per patch; (c) median inters. area b/w clusters of lymp. and non-lymp.; (d) side length min/max (Delaunay triang.) of lymp. clust.; (e) skewness of compactness of lymphocyte clusters; and (f) triangle area min/max (Delaunay) of lymphocyte clusters.

3.2.

Experiment 2: Comparing Patterns of Immune Response in COVID-19 and H1N1

A ROC-based feature selection algorithm⁴² was employed to find the 10 most discriminating features. The Wilcoxon test was employed to determine whether differences with respect to the top features were statistically significant. This process was also carried out for four features related to lymphocytes, namely quantity of lymphocytes, area of lymphocyte clusters, and compactness of cluster centroids. Since multiple comparisons were made at the patch level and there were many features per patient, the $p$ values corresponding to a comparison of image patches between COVID-19 versus H1N1 were adjusted using the Benjamini-Hochberg procedure using a false discovery rate (FDR) of 10% (Table 1). As all the assessed $p$ values were smaller than the critical value, they were considered significant.

Table 1

Benjamin-Hochberg procedure using a FDR of 10% for the top discriminating features (as dictated by an ROC-based feature selection algorithm) and four lymphocyte-related features (quantity of lymphocytes, area of lymphocyte clusters, and compactness of cluster centroids).

Variable	p value	Rank	Adjusted p value using FDR approach
Top feat 1	$2.12\times {10}^{-25}$	1	0.0071
Top feat 2	$7.24\times {10}^{-19}$	2	0.0143
Top feat 5	$2.65\times {10}^{-18}$	3	0.0214
Top feat 6	$2.94\times {10}^{-18}$	4	0.0286
Top feat 3	$4.58\times {10}^{-16}$	5	0.0357
Top feat 4	$6.27\times {10}^{-16}$	6	0.0429
Top feat 7	$4.48\times {10}^{-15}$	7	0.0500
Top feat 9	$6.16\times {10}^{-15}$	8	0.0571
Top feat 10	$2.19\times {10}^{-14}$	9	0.0643
Top feat 8	$1.82\times {10}^{-13}$	10	0.0714
Lymp feat 4	$1.84\times {10}^{-5}$	11	0.0786
Lymp feat 3	$5.02\times {10}^{-5}$	12	0.0857
Lymp feat 2	$1.44\times {10}^{-4}$	13	0.0929
Lymp feat 1	$2.22\times {10}^{-4}$	14	0.1000

Figure 6 illustrates the top four most discriminating features, which are related to overlap between clusters of lymphocytes and non-lymphocytes. Although interpretation of these features is subtle, they indicate that very large sized clusters of non-lymphocytes are more frequently found in H1N1 compared to COVID-19. Figure 7 illustrates the four lymphocyte-related features analyzed. We can deduce that in COVID-19 there are more clusters of lymphocytes, of smaller size, that are closer to each other (high compactness) compared to H1N1, lymphocyte clusters tending to be fewer and larger.

Fig. 6

Split violin plots comparing the variation of the 4 top discriminating features in patches of COVID-19 (9 patients, 496 patches) and H1N1 (2 patients, 110 patches). These features are related to the intersection of clusters of lymphocytes and non-lymphocytes and are defined as follows. (a) $\frac{AF1\cap AF2}{AF2}$, where $AF1$ is the area of a convex hull containing the centroids of all the lymphocyte clusters and $AF2$ is the area of a convex hull containing the centroids of all the non-lymphocyte clusters. (b) Total $2*\frac{A1\cap A2}{A1+A2}$, where $A1$ is the area of clusters of lymphocytes and $A2$ is the area of clusters of non-lymphocytes. (c) $2*\frac{AF1\cap AF2}{AF1+AF2}$. (d) $AF1\cap AF2$.

Fig. 7

Split violin plots comparing the variation of lymphocyte-related features in patches of COVID-19 (9 patients, 496 patches) and H1N1 (2 patients, 110 patches) patients. (a) Median compactness of lymphocyte clusters per patch; (b) number of clusters of lymphocytes per patch; (c) number of lymphocytes per patch; and (d) max area of lymphocyte clusters per patch.

Figure 8 shows spatial heat maps illustrating the variation of the number of non-lymphocytes per patch across 2 WSIs, one for a patient who died of COVID-19 (OCME) and another of H1N1 (CC). The density of non-lymphocytes per patch is higher for H1N1.

Fig. 8

Spatial heat maps illustrating the number of non-lymphocytes per patch and representative patches of a patient who died of COVID-19 and another who died of H1N1. $H$ stands for a high number of non-lymphocytes and $L$ for low. The patches in the bottom row show a hot spot for each disease (i.e., containing a high number of non-lymphocytes), revealing that in H1N1 there is a higher density of non-lymphocyte nuclei (green dots) compared to COVID-19.

Finally, t-SNE was run to evaluate whether spatial-architecture-related features were able to separate H1N1 and COVID-19 patches. Figure 9(a) shows no evident clusters across the different sites, suggesting that data are similarly distributed. Figure 9(b) presents the distribution of patches according to diagnosis, exhibiting an area with a high concentration (top middle) of H1N1 points and a subtle grouping of H1N1 points (left). This could imply that some architectural patterns are more commonly found in H1N1 than in COVID-19. When looking at the architectural patterns of the studied nuclei families, in general, lymphocyte clusters are larger in H1N1. Similarly, non-lymphocytes are more numerous and compact, tending to generate larger clusters than COVID-19 (Fig. 9, red/yellow-border patches).

Fig. 9

Unsupervised clustering of spatial architectural features of immune clusters of COVID-19 and H1N1 patches using t-SNE and illustration of representative patches. (a) Each point represents an image patch, and each color corresponds to a different institution. The t-SNE algorithm started out with 350 features and reduced these to 2. (b) The same clustering but in this case each point represents a patch and each color a disease (i.e., either H1N1 or COVID-19). The red-border patches correspond to H1N1 and exhibit large clusters of non-lymphocytes (green border), which suggest that these nuclei are numerous and compact. The yellow-border patches correspond to COVID-19 and present as smaller groups of non-lymphocytes.

4.

Discussion and Conclusion

COVID-19 causes a wide spectrum of signs and symptoms.^43,44 In severely ill cases, some authors mentioned an exuberant cytokine production along with an intrinsic viral capacity to limit the innate immune response, which causes severe organ damage and perpetuates viral replication.^{1,2,19,45–48} Previous works have employed peripheral blood for describing the immune behavior in COVID-19 patients. The most common findings include low counts of total number of T lymphocytes and their CD4 subtype.^{15,46,49–51} These lymphocyte counts have been even proposed as predictors of outcome.^49–51

In tissue, visual pathologic assessment of lymphocytes in autopsied samples of patients who have died of COVID-19 has been reported. This estimation of lymphocytes is, however, time-consuming^26,52 and has been subjectively estimated as: low-versus-medium, mild-versus-moderate, or patchy-versus-diffuse accompanying changes such as diffuse alveolar damage, acute bronchopneumonia, and microthrombi.^6,7,10,12,53 No comparison of immune cell response among different stages in lung tissue of COVID-19 or other viral diseases has been performed to date.

To our knowledge, this is the first study in: (1) describing quantitatively immune nuclei patterns in lung tissue from patients who died of COVID-19 and (2) exploring the differences in the morphologic immune architecture of patients who succumbed to COVID-19 and H1N1, a comparable viral infection whose clinical spectrum, laboratory test results, pharmacotherapy, and prevention strategies have been extensively explored.^15,54

Using an ROC-based feature selection method, violin plots, and t-SNE on patches of COVID-19 and H1N1, we found that spatial arrangement of lymphocytes per tile was different in both diseases (Figs. 5, 6, and 9). In COVID-19, we found three important patterns: (1) multiple smaller-sized lymphocyte clusters, (2) immune clusters typically surrounded by other lymphocyte clusters, and (3) regions encompassed by the lymphocyte clusters typically overlapping with other non-lymphocyte clusters. These patterns appear to suggest that the lymphocyte response in patients who died of COVID-19 is relatively modest. A possible explanation for this phenomenon is the widely reported lymphocytopenia;^49–51 however, a more focused study needs to be performed to assess any correlation between lymphocytic infiltration in lungs and lymphocyte levels in the blood.

In H1N1 patients, by contrast, we noticed (1) fewer and larger clusters of lymphocytes and (2) significant presence of few but extremely large clusters of non-lymphocytes. This appears to be concordant with previous findings reporting a considerable quantity of macrophages accompanying the lymphocytic infiltration in H1N1, as well as a higher prevalence of acute immune cell infiltration.^21,22 By contrast, fewer neutrophils have been found in COVID-19 cases given the potential of the $N$ protein (found in SARS-CoV-2) to block IFN-1 production, a cytokine involved in neutrophil attraction.^45,55 This imbalanced proportion of immune cells in solid tissues such as the lung is a potential explanation for the differences we found when comparing lymphocyte spatial distribution and their interaction with other nuclei in COVID-19 versus H1N1 patients.

In fatal cases of COVID-19, the inability of lymphocytes to limit viral damage despite possibly being extremely activated could be a consequence of abnormal lymphocyte–lymphocyte or lymphocyte–non-lymphocyte crosstalk, concordant with a prior hypothesis by Tan et al.⁵¹ In contrast, in H1N1, the response of non-lymphocytes appears to be flourishing, possibly reflecting either better coordination between lymphocytes and their surrounding nuclei, or a key role of superimposed acute inflammation.

We acknowledge the limitations of this work. Quantitative differences in a small number of cases might not hold up when larger numbers of cases are studied; however, it is worth mentioning that studies involving autopsies usually draw conclusions from fewer cases. Calabrese et al.⁵⁶ presented a work summarizing the pathological findings described in 23 articles reporting autopsies in suspected/known COVID-19 patients, evidencing that the median number of cases per study is 2 patients. Of course, as more cases are published, a more comprehensive analysis of the spatial architecture of the immune response for COVID-19 and H1N1 could be performed. Some studies⁶ have employed immunohistochemistry and have been able to analyze the importance of subfamilies of lymphocytes (e.g., CD4 or CD8); however, in those works, the spatial architecture of such families was not analyzed. Our findings using standard H&E samples suggest that spatial relationships between immune and non-immune nuclei could play a role in characterizing these diseases. Identification of subgroups of lymphocytes or other immune nuclei could meaningfully enrich a similar work in the future. Future work may also benefit from including an extensive comparison against effective immune response as well as with normal lymphoid tissue associated to airway mucosa. Additionally, a comprehensive analysis of the implications of the threshold value employed to build nuclei clusters for characterizing both diseases could meaningfully enrich this work in the future.

In summary, we used a computational-pathology-based image analysis approach to characterize the immune nuclei architecture in digitized autopsied H&E lung samples from patients who died of COVID-19 and H1N1. As far as we know, this is the first attempt to quantitatively capture morphology for characterizing lymphocytic infiltrates and their relationship to surrounding nuclei in COVID-19 and simultaneously compare this arrangement in H1N1. In this preliminary study, we found that the immune response in COVID-19 consistently showed multiple lymphocyte clusters of small size, suggesting that the lymphocyte response is rather modest, possibly due to lymphocytopenia. In contrast, in H1N1, we found larger lymphocyte clusters accompanied by large clusters of non-lymphocytes, a possible reflection of increased prevalence of macrophages, and other immune cells. In combination with cytokine profiling⁵⁷ and after testing this algorithm in a larger population, our approach could potentially help enrich and inform understanding in the process of developing therapeutic involving immune response modulation in COVID-19 patients.