Early assessment of lung function in coronavirus patients using invariant markers from chest X-rays images

The primary goal of this manuscript is to develop a computer assisted diagnostic (CAD) system to assess pulmonary function and risk of mortality in patients with coronavirus disease 2019 (COVID-19). The CAD system processes chest X-ray data and provides accurate, objective imaging markers to assist in the determination of patients with a higher risk of death and thus are more likely to require mechanical ventilation and/or more intensive clinical care.To obtain an accurate stochastic model that has the ability to detect the severity of lung infection, we develop a second-order Markov-Gibbs random field (MGRF) invariant under rigid transformation (translation or rotation of the image) as well as scale (i.e., pixel size). The parameters of the MGRF model are learned automatically, given a training set of X-ray images with affected lung regions labeled. An X-ray input to the system undergoes pre-processing to correct for non-uniformity of illumination and to delimit the boundary of the lung, using either a fully-automated segmentation routine or manual delineation provided by the radiologist, prior to the diagnosis. The steps of the proposed methodology are: (i) estimate the Gibbs energy at several different radii to describe the inhomogeneity in lung infection; (ii) compute the cumulative distribution function (CDF) as a new representation to describe the local inhomogeneity in the infected region of lung; and (iii) input the CDFs to a new neural network-based fusion system to determine whether the severity of lung infection is low or high. This approach is tested on 200 clinical X-rays from 200 COVID-19 positive patients, 100 of whom died and 100 who recovered using multiple training/testing processes including leave-one-subject-out (LOSO), tenfold, fourfold, and twofold cross-validation tests. The Gibbs energy for lung pathology was estimated at three concentric rings of increasing radii. The accuracy and Dice similarity coefficient (DSC) of the system steadily improved as the radius increased. The overall CAD system combined the estimated Gibbs energy information from all radii and achieved a sensitivity, specificity, accuracy, and DSC of 100%, 97% ± 3%, 98% ± 2%, and 98% ± 2%, respectively, by twofold cross validation. Alternative classification algorithms, including support vector machine, random forest, naive Bayes classifier, K-nearest neighbors, and decision trees all produced inferior results compared to the proposed neural network used in this CAD system. The experiments demonstrate the feasibility of the proposed system as a novel tool to objectively assess disease severity and predict mortality in COVID-19 patients. The proposed tool can assist physicians to determine which patients might require more intensive clinical care, such a mechanical respiratory support.

, and data from the University of Louisville, USA and Mansoura University, Egypt. The research protocol was approved by the institutional review board (IRB) at the University of Louisville and Mansoura university as well as all methods were performed in accordance with the relevant guidelines and regulations and the patients informed consent was obtained. For the patients who passed away, an informed consent was obtained from legal guardian/Next of kin for the dead patients. These databases include 200 subjects tested as COVID-19 positive, 100 from patients who eventually died from the infection and 100 patients who ultimately recovered. These databases comprise a heterogeneous collection of digital X-ray images, which was the primary motivation to develop rotation, scale, and translation invariant MGRF model from which we extract the proposed imaging markers to grade the severity of lung infection in COVID-19 patients. We did our experiments using the merged datasets from all three databases to overcome the issue of data balance as the number of subjects for the two classes, (high and low severity cases), in each database was different. The dead cases had been confirmed based on the following radiology protocol 33 : (i) Bilateral and predominantly peripheral opacifications and/or consolidations were rated as typical for a COVID-19 infection. (ii) A distribution pattern with opacifications and/or consolidations limited to one pulmonary lobe consistent with lobar pneumonia was rated as non-typical for a COVID-19 infection.
All changes that could not be classified as non-typical or typical were rated as indeterminate. In a subgroup of these patients with indeterminate findings, soft criteria for a possible COVID-19 infection were defined as the unilateral presence of predominantly peripheral.
Proposed computer aided diagnostic (CAD) system. The proposed CAD system to detect the severity of lung infection is shown in Fig. 1. The CAD system consists of three major steps: (i) preprocessing steps to improve contrast of the X-ray images and identify the region of interest in order to enhance diagnostic accuracy of subsequent steps; (ii) modeling the appearance of infected chest tissue using a new Markov-Gibbs random field (MGRF) constructed to be invariant under rotation, translation, and change of scale; and (iii) a neural network (NN)-based fusion and diagnostic system to determine whether the grade of lung infection is low or high.
Data preprocessing. To improve the accuracy of the proposed approach, we manually segmented the lung region from the original X-ray image, Fig. 2a, as demonstrated in Fig. 2b. The second step is to enhance lung tissue contrast, for which we use regional dynamic histogram equalization (RDHE) 34 . Proper analysis of the type of noise present in the chest X-ray image may help to select proper denoising methods, which preserve the important texture information while reducing the noise 35,36 . The RDHE approach divides the image into blocks x rows high by y columns wide. Then, dynamic histogram equalization is applied within each block to adaptively enhance the contrast. Therefore, the image histogram is remapped by block, and pixel values are adjusted relative to the other pixels in their x × y neighborhood. The contrast-enhanced X-ray image resulting from the RDHE approach is illustrated in Fig. 2c.  www.nature.com/scientificreports/ The third preprocessing step is to identify and mask off the healthy lung tissues from the infected tissues. This step narrows the search space to focus only on the abnormal tissues and serves to increase the diagnostic accuracy of the CAD system. To achieve this step, we use our previously published methodology 37 that considers both the spatial interaction between nearby image pixels and the intensity distribution of those pixels within the lung region of interest. We follow the conventional description of the MGRF model in terms of independent signals (images) and interdependent region labels (segmentations); yet we focus on more accurate model identification 37 . Each image segment corresponds to a single dominant mode of the empirical distribution (i.e. histogram) of gray levels. To identify the dominant modes, each image histogram is considered to be sampled from a linear combination of discrete Gaussians (LCDG) distribution 37 . We fit an initial LCDG model to the empirical distribution using a modified expectation-maximization (EM) algorithm 37 . Free parameters of the LCDG to be optimized are the number of discrete Gaussian components and their respective weights (positive and negative), shape, and scale parameters. Then, the initial LCDG-based segmentation is iteratively refined using the MGRF model with its analytically estimated potentials 37 . Figure 2d shows the extracted pathological tissues using our proposed algorithm. Additional details can be found in El-Baz et al. 37 .
Rotation, scale, and translation invariant MGRF model. We constructed the proposed MGRF model in order that the image need not be aligned with any particular frame of reference in order to use it to grade the severity of lung infection (low vs. high). To construct the appearance of the infected lung regions, we consider the X-ray images as samples from a piecewise stationary MGRF with a central-symmetric system of pixel-pixel interactions. Let n ν denote a set of central-symmetric pixel neighborhoods indexed by ν ∈ {1, . . . , N} . Each n ν is a set of coordinate offsets (ξ , η) specified by a semi-open interval of interpixel distances (d ν,min , d ν,max ] such that the n ν -neighborhood of pixel (x, y) comprises all pixels (x ′ , y ′ ) such that d ν,min < (x − x ′ ) 2 + (y − y ′ ) 2 ≤ d ν,max . A neighborhood system corresponding to d ν,min = ν − 1 2 and d ν,max = ν + 1 2 , ν ∈ {1, 2, 3} , is illustrated in Fig. 3. Associated with the neighborhood system is a set of N + 1 Gibbs potential functions of gray value and gray value co-occurrences V 0 : Q → R and V ν : Q × Q → R , ν ∈ {1, . . . , N} , where Q is the range of pixel gray levels, e.g. Q = {0, . . . , 255} in the case of 8-bit images.
For a given image/label map pair (g t , m t ) from our training set S, t ∈ {1, . . . , T} , let R t = {(x, y) | m t (x, y) = ob} denote the subset of the pixel lattice supporting the infected lung region. Denote the set of n ν -neighboring pixels restricted to R t by  www.nature.com/scientificreports/ Finally let f 0,t and f ν,t , ν ∈ {1, . . . , N} denote empirical probability distributions (i.e., relative frequency histograms) of gray values and gray value co-occurrences in the training infected region from the X-ray image g t , The joint probability of object pixels in image g t according to the MGRF model is given by the Gibbs distribution where ρ ν,t = |C ν,t |/|R t | is an average cardinality of n ν over the sublattice R t .
Assuming lungs having the same pathology exhibit similar morphology in X-ray images, then we may approximate the previous expressions by their averages over the training set S: |R t | ≈ R ob and |C ν,t | ≈ C ν,ob . Here If we assume further that the observations in S are statistically independent (e.g., each X-ray is taken from a different patient), the expression for joint probability of object pixels may be likewise simplified 38 : Here, ρ ν = C ν,ob /R ob , and the probability vectors F pix,ob and F ν,ob are the averages of the relative frequency histograms and normalized gray level co-occurrence matrices, respectively, over all objects in the training set. The problem of zero empirical probabilities, which can arise when a relatively small volume of the training data is available to identify the MGRF model, is dealt with using pseudocounts. Then Eqs. 1 and 2 are modified as follows The Bayesian quadratic loss estimate suggests using the offset ε = 1 , while a more conservative approach 38 suggests using ε = 1/Q in Eq. 4 and ε = 1/Q 2 in Eq. 5.
Using the same analytical approach as in Ref. 38 , the Gibbs potential functions are approximated using the centered, training-set average, normalized histograms and co-occurrence matrices: Using the above estimated potentials, we can calculate the Gibbs energy of the infected lung region b in an X-ray image g as follows: Here, N ′ is a selected top-rank index subset of the neighborhoods, and the empirical probability distributions F 0 and F ν are calculated over the object pixels b of g.
To summarize, the whole training approach is as follows: 1. Read all infected regions from the training data having class "severe" lung infection.

NN-based fusion and diagnostic system
A new NN system that can fuse the diagnostic results from the three estimated Gibbs energy at three different radii is developed. The proposed NN-based model consists of four blocks as illustrated in Fig. 4. Three of them are fed with the three different cumulative distribution functions (CDFs) of the estimated Gibbs energy, then (1) (3) In order to tune the hyper-parameters used in our proposed NN system, a hyper-parameters estimation approach is used. The parameters to be estimated are the number of bins used to calculate CDF, number of hidden layers in the NN model, number of neurons in the hidden layer, and finally the activation function used to calculate the output of each neuron. We ran several experiments using random values for these parameters to estimate their optimal values using training data. All the results that are demonstrated in the "Experimental results" section have been obtained using the following setting: to handle all energy values, the chosen value for the number of CDF bins is 175; the number of hidden layers in NN 1, NN 2, and fusion NN is one, while for NN 3, there are no hidden layers (searching from 0 to 10); the number of neurons per hidden layer is 50, 20, and 2 for NN 1, NN 2, and Fusion NN, respectively (searching from 1 to 200); and finally, the sigmoid activation function has been selected after also considering the tangent and softmax activation functions.

Experimental results
To highlight the innovation in our approach, we demonstrate the Gibbs energy calculated at three radii as a color map fused over the X-ray images. One example for which it holds, it is clear from Fig. 5 that the Gibbs energy in cases of high severity of COVID-19 pneumonia is high compared with the Gibbs energy for low-severity COVID-19 pneumonia. Since the collected X-ray images have different resolutions, we use CDF as a new scale-invariant representation to the estimated Gibbs energy which makes it suitable for all data collection protocols as shown in Fig. 6. To highlight the advantage of the proposed Gibbs energy as a new discriminatory image marker, we calculate the average CDFs with a demonstration of the standard deviation at each point for both classes (high severity vs. low severity). As is clear from Fig. 7, the CDFs are rather distinctive which allows for straightforward classification by the proposed NN-based classifier. The output of the CAD system was an assessment of the severity of pneumonia in COVID-19 patients with two levels: a low severity of infection ("low") or high severity of infection ("high") as shown in Fig. 4. This was compared to the ground truth of the 200 clinical cases collected, 100 of which were from patients who died of COVID-19 and 100 of which recovered. Accurate system outputs include an assessment of "low" in a case that recovered and an output of "high" in a case that died. To confirm the accuracy of the proposed NN classification and fusion system, leave-one-subject-out (LOSO), tenfold, fourfold,  Table 1. We use the following objective metrics to measure the accuracy of the proposed NN-based fusion system: (i) sensitivity, (ii) specificity, (iii) accuracy, and (iv) Dice similarity coefficient (DSC). As demonstrated in Table 1, the proposed system has achieved 100% accuracy with the LOSO validation test and 98.00% ± 2.00% for a twofold validation test (real-life scenario), all of which confirm the high accuracy of the proposed CAD system. To highlight the contribution of each Gibbs energy at each radius, we construct an NN-based classifier using the estimated Gibbs energy at each radius. As is clear from Table 1, the NN-classifier based on the estimated   www.nature.com/scientificreports/ Gibbs energy at ν 3 demonstrates the highest accuracy compared with the classification accuracies based on the estimated Gibbs energy at ν 2 and ν 1. Also, it is worth mentioning that fusing the three estimated Gibbs energies by using the NN-Based classification system achieves higher accuracy compared with classification accuracies based on each single estimated Gibbs energy. Finally, to highlight the accuracy of the proposed NN-based fusion system, we compare its accuracy with support vector machine (SVM), random forest, naive Bayes, K-nearest neighbors (KNN), and decision trees classifiers. Table 2 clearly shows that the NN-based classification and fusion system has achieved the highest accuracy compared with other approaches. Statistical significance of the choice of feature set or classifier architecture on diagnostic performance was assessed using Friedman rank sum tests [39][40][41][42] . Post hoc comparison of our proposed NN-based classifier with the alternatives was done using Wilcoxon signed rank tests with Bonferroni correction to the estimated p−values. Feature selection had a significant impact on classifier performance (Table 1), with Friedman test χ 2 = 32.7 , 3 d.f., p = 3 × 10 −7 . The fused feature set was shown to be preferable to v1 (Wilcoxon V = 118 , Bonferroni corrected p = 0.0033 ), v2 ( V = 136 , p = 0.0014 ), and v3 ( V = 136 , p = 0.0014 ). Different classifier architectures also produced significantly different results (Table 2)

Discussion and conclusion
ARDS is the most common and severe pulmonary complication in COVID-19 patients. It is an acute hypoxemic respiratory failure that requires oxygen and ventilation therapy including intubation and invasive ventilation. Clinically patients may have dyspnea, tachypnea (respiratory rate ≥ 30 breaths per minute), decreased peripheral oxygen saturation SpO 2 ≤ 93% , poor oxygenation with the ratio of the partial pressure of arterial oxygen to fraction of inspired oxygen PaO 2 /FiO 2 < 300 mmHg, or lung infiltrates greater than 50% within 48 h 9 . ARDS occurred in 20% of hospitalized patients and 61% of ICU patients in Zhongnan Hospital in Wuhan 3,4 . ARDS occurs when capillaries in the lung leak fluid into the alveoli, thereby impairing gas exchange in the lung and reducing oxygen uptake into the systemic arterial circulation. The consequent decrease in blood oxygen levels can be directly life-threatening, leading to multi-organ failure. Respiratory support of COVID-19 may use invasive or non-invasive methods to force oxygen into the airways under pressure. Invasive ventilation uses an endotracheal tube to feed oxygen directly into the lungs. Non-invasive methods employ such devices as continuous positive airway pressure (CPAP) and oxygen hoods; there is no use of an internal tube, and they are used in the management of less severe cases.
Despite being vital for supporting respiration in patients with ARDS, ventilators are in short supply in hospitals. According to Imperial College London, 30% of patients diagnosed with COVID-19 are strongly recommended to be admitted to the hospitals, with a significant fraction of those patients also requiring respiratory support. As the pandemic spreads across the world, many countries stopped exporting ventilators 43,44 . The paucity of ventilators is even more acute in under developed and developing countries in South America, Asia, and Africa. www.nature.com/scientificreports/ High-pressure ventilation may cause lung injury, also called barotrauma or ventilator-induced lung injury (VILI). Even non-invasive ventilation carries some risk, as stress and strain associated with high tidal volumes may cause patient self-induced lung injury (P-SILI). The additional inflammation due to VILI or P-SILI may lead to aggravation of pulmonary edema and worsening of the very respiratory distress that ventilation was intended to treat. There is also the risk of heart failure, hypervolemia, and multi organ dysfunction, alone or in combination 45 . Unfortunately, COVID-19 patients who are admitted to the ICU and require mechanical ventilation show strikingly high rates of mortality, ranging from 50 to 97% early in the pandemic [46][47][48][49][50][51] . A more recent study from Emory University showed lower but still dramatic mortality rates of 36% in ICU patients requiring mechanical ventilation and 30% in all COVID-19 patients admitted to the ICU 52 .
Accurate and rapid diagnosis of COVID-19 pneumonia severity is very challenging for radiologists as the disease has rapidly spread across the globe. Based on the results demonstrated in this study, AI systems, especially those based on deep learning, are promising tools to assist initial screening by radiologists. It could decrease workload, improve diagnostic accuracy, and enable appropriate treatments and ventilation management of COVID-19 patients. In the case of a pandemic as we now face, medical resources are seriously strained and must be used as efficiently as possible. Rapid diagnosis and accurate prognosis are essential. The AI-based method shows great potential to quantify disease severity and could be used to inform treatment decision-making in patients with COVID-19. AI in concert with thoracic imaging and other clinical information (epidemiology, PCR, clinical symptoms, and laboratory indicators) can effectively improve clinical outcomes 53 . AI can increase the utility of chest X-ray imaging beyond first-line diagnostic imaging and into the areas of risk stratification, monitoring of clinical course, and selection between management approaches, such as invasive vs. non-invasive ventilation, for COVID-19 patients. Multimodal data, be they clinical, epidemiological, or potentially molecular data, can by fused with imaging data in an AI framework to build systems to detect and treat COVID-19 patients and potentially to contain its spread 54 . Moreover, we plan to work on X-ray scans/data that are collected at different time points to evaluate the progressing of the infection/pneumonia with the treatment course.
In conclusion, the results herein demonstrate the feasibility of using AI with chest X-ray imaging data to determine the severity of lung involvement in cases of COVID-19. Severity of pneumonia on chest X-ray correlated highly with mortality in this study, and thus this CAD system can potentially also be used to predict mortality in COVD-19 patients.

Data availability
Materials, data, and associated protocols will be available to readers after the manuscript being accepted.