Predicting multiple sclerosis disease severity with multimodal deep
  neural networks
Multiple Sclerosis (MS) is a chronic disease developed in human brain and
spinal cord, which can cause permanent damage or deterioration of the nerves.
The severity of MS disease is monitored by the Expanded Disability Status Scale
(EDSS), composed of several functional sub-scores. Early and accurate
classification of MS disease severity is critical for slowing down or
preventing disease progression via applying early therapeutic intervention
strategies. Recent advances in deep learning and the wide use of Electronic
Health Records (EHR) creates opportunities to apply data-driven and predictive
modeling tools for this goal. Previous studies focusing on using single-modal
machine learning and deep learning algorithms were limited in terms of
prediction accuracy due to the data insufficiency or model simplicity. In this
paper, we proposed an idea of using patients' multimodal longitudinal and
longitudinal EHR data to predict multiple sclerosis disease severity at the
hospital visit. This work has two important contributions. First, we describe a
pilot effort to leverage structured EHR data, neuroimaging data and clinical
notes to build a multi-modal deep learning framework to predict patient's MS
disease severity. The proposed pipeline demonstrates up to 25% increase in
terms of the area under the Area Under the Receiver Operating Characteristic
curve (AUROC) compared to models using single-modal data. Second, the study
also provides insights regarding the amount useful signal embedded in each data
modality with respect to MS disease prediction, which may improve data
collection processes.
Kai Zhang, John A. Lincoln, Xiaoqian Jiang, Elmer V. Bernstam, Shayan Shams
# Predicting multiple sclerosis disease severity with multimodal deep neural networks

###### Abstract

Multiple Sclerosis (MS) is a chronic disease developed in human brain and spinal cord, which can cause permanent damage or deterioration of the nerves. The severity of MS disease is monitored by the Expanded Disability Status Scale (EDSS), composed of several functional sub-scores. Early and accurate classification of MS disease severity is critical for slowing down or preventing disease progression via applying early therapeutic intervention strategies. Recent advances in deep learning and the wide use of Electronic Health Records (EHR) creates opportunities to apply data-driven and predictive modeling tools for this goal. Previous studies focusing on using single-modal machine learning and deep learning algorithms were limited in terms of prediction accuracy due to the data insufficiency or model simplicity. In this paper, we proposed an idea of using patients' multimodal longitudinal and longitudinal EHR data to predict multiple sclerosis disease severity at the hospital visit. This work has two important contributions. First, we describe a pilot effort to leverage structured EHR data, neuroimaging data and clinical notes to build a multi-modal deep learning framework to predict patient's MS disease severity. The proposed pipeline demonstrates up to 25% increase in terms of the area under the Area Under the Receiver Operating Characteristic curve (AUROC) compared to models using single-modal data. Second, the study also provides insights regarding the amount useful signal embedded in each data modality with respect to MS disease prediction, which may improve data collection processes.

Database; Deep neural network; Multiple sclerosis; Expanded disability status scale

## Background

Recent advantages in deep learning have shown success in various areas of healthcare, such as brain Magnetic Resonance Imaging (MRI) automatic volume segmentation and classification [1], clinical text mining and disease prediction [2], risk predictions [3], etc. The fast-growing Electronic Health Records (EHR) in healthcare provides a great number of opportunities for both the data mining and deep learning communities to explore the rich information embedded in different data modalities and tap the potentiality of using this information for predictive modeling, to benefit effective healthcare delivery and better-quality caring for patients.

Multiple sclerosis (MS) is a potentially disabling disease that affects the human brain and spinal cord. An estimation of MS prevalence by the year 2010 of 10-year accumulation shows there are over 700,000 MS cases in adults in the United States [4]. Recent advantages in MS disease research found that patients who died from MS suffer up to 39% neuron count loss compared to usual patients without MS [5]. The
human brain has mechanisms for self-repair and regenerative potential that could repair the brain plaques [6], however, such ability is very limited. Therefore, prompt action to prevent or slow down brain damage is critical to MS disease treatment [7]. Effective treatment relies on a correct grading of the MS severity, and scoring systems are widely used to achieve this goal. The Expanded Disability Status Scale (EDSS) score [8] is a widely used ordinal scoring system by healthcare providers to monitor clinical disability in MS. It is composed of diverse functional systems, including pyramidal functions (muscle strength, tone, and reflexes), cerebellar functions (coordination and balance), brainstem functions (eye movements, speech, and swallowing), sensory functions (light touch, pain, and vibratory sense), bowel and bladder functions, visual functions, cerebral (cognition) and ambulation. Based on EDSS, Roxburgh et al. proposed a Multiple Sclerosis Severity Score (MSSS) which can be used to determine MS disease progression using single assessment data (when a patient has only one assessment during the disease course) [9].

Several milestones of the EDSS score have been commonly used to define different stages of the MS disease course. The EDSS 4 (significant disability but able to walk without aid or rest for 500 m), EDSS 6 (requires unilateral assistance to walk about 100 m with or without resting) and EDSS 7 (ability to walk no more than 10 m without rest while leaning against a wall or holding onto furniture for support) were commonly-used milestones for studying MS disease severity. For example, Confavreux et al. used the above milestones to study the effect of relapses on the progression of irreversible disability [10]. The same milestones have also been used to study the contribution of relapses to worsening disability and evaluate the MS therapies' effect on delaying the disability accumulation [11]. A Sweden research group studied whether the risk of reaching the above disability milestones in MS has changed over the last decade [12]. Rzepinski et al. used the EDSS milestones to explore early clinical features of MS and how they affect patients' long-term disability progression [13]. The same milestones were also used to study how these factors affect the time to transition from relapsing-remitting MS (RRMS) to secondary progressive MS (SPMS).

A patient's EDSS score needs to be evaluated by a well-trained specialist to ensure that the assessment is correctly performed, which limits the application of EDSS to clinics with MS disease specialties. Several research studies have attempted to address this problem using machine learning or deep learning models. In particular, Pinto et al. proposed to use machine learning models to predict MS progression, based on the clinical characteristics of the first five years of the disease [14]. Zhao et al. used a support vector machine (SVM) classifier and demographic, clinical, and MRI data obtained at years one and two to predict patients' EDSS at five years follow-ups [15]. Sacca et al. explored different machine learning models (Random Forest, Support Vector Machine, Naive-Bayes, K-nearest-neighbor, and Artificial Neural Network) and used the features extracted from functional-MRI to perform MS disease severity classification [16]. Narayana et al. proposed to use the VGG-16 convolutional neural network (CNN) to predict enhancing lesions in MS patients using non-contrast MRIs [17]. D'Costa et al. proposed a transformer model named MS-BERT to predict EDSS score from patient's neurological consult note [2]. Ciotti proposed a clinical instrument to retrospectively capture levels of EDSS and the algorithm got a Kappa score of 0.80 between captured EDSS and real EDSS [18].
Chase et al. also used neurological consult notes but with simpler models (Naive Bayes classification model) and features (word frequency) [19]. Dekker et al. used multiple linear regression models on patient brain lesion volumes and its variation over the years to predict physical disability [20]. The aforementioned studies explored the idea of using machine learning and deep learning methods on various modalities of EHR datasets to predict patient's EDSS of the current hospital visit or in the near future. The above works only explored a limited amount of patient information (either clinical notes, or basic lesion volume information extracted from MRI, or patient clinical characteristics), by adopting the off-the-shelf machine learning models, or deep learning models that were developed for general tasks and utilized without being customized to capture the complex nature of this prediction problem. Based on the above studies and the recent research advances in multimodal deep learning, it is reasonable to assume that using multimodal deep learning methods could integrate fragmented information from each modality and brings more accurate predictions for MS disease. Therefore, this study tried to answer the question of whether can we harmonize all the available EHR data modalities collected from patient clinic visits and use longitudinal data to perform more accurate MS severity prediction. A few research study findings have found that MRI data and some lab tests can contain useful information about MS disease severity. For example, studies have shown that the thickness of cortical and deep grey matter has a high correlation with the MS disease severity, suggesting that the MRI images are an informative data source to predict MS severity [21, 22]. Some laboratory tests were also documented as playing an important role in this regard, such as the cerebrospinal fluid (CSF) [23, 24], and serum neurofilament light chain (nFl) [25].

This study tried to answer the above question using a data-driven approach. We explored the idea of using patients' MRI images, clinical notes, and structured EHR data (including laboratory tests, vital sign observations, medication prescriptions, patient demographics) that were collected during patients' clinic visits to predict MS disease severity at the visit. We propose a multimodal deep neural network that takes MS patient Electronic Health Records of multiple modalities, including the MRI images (pre- and post-contrast T1 weighted image, T2 weighted image, fluid-attenuated inversion recovery image and proton density image), patient's clinical notes data and structured EHR (laboratory tests, vital signs, medications, demographic information) to perform MS disease severity prediction.

We also propose to use patients' longitudinal data for EDSS milestone prediction, based on the fact that evidence about patients' MS disease severity should not only be embedded in the most recent EHR data but also richly contained in data of all previous clinic visits. Compared to using cross-sectional data (e.g. using clinical notes of the current visit to predict EDSS score [2]), we propose to use the patient's both current clinic visit and historical EHR data to train a multimodal deep neural network. The longitudinal data contains more MS disease progression information compared to cross-sectional data and will help the model make more accurate predictions of the patient's status at the moment. The contributions of this study are three-fold.
* A novel deep learning architecture (a multimodal neural network) and data fusion mechanism which takes Electronic Health Records including medications, vital signs, laboratory test results, clinical imaging, and physician notes to tackle the difficult problem of MS disease severity prediction. The results show significant prediction accuracy improvements compared to using single-modality data or simpler models.
* Using longitudinal data (both current and historical visits data) instead of cross-sectional data (data of current visit) to accurately classify patient EDSS score milestones at the current clinic visit.
* Exhibits how much useful information is embedded in each data modality for the prediction of MS severity. Various attention mechanisms are adopted in the proposed neural network to provide model explainability and enhance prediction accuracy.
* An end-to-end AI model that works on readily available data with a limited pre-processing process (e.g. does not need feature extraction as a pre-processing step, such as extracting the thalamic volume, lateral ventricle volume, etc. to train the model)

The paper is structured as follows. The Data section explains the dataset that was used in this study and the details of each data modality. The subsequent section explains our designed deep neural network architecture, followed by an Experiment section including our experiment design, the obtained results, and the discussion. Finally, a summary and conclusion are given in the Discussion section.

### Data

Our database contains a rich set of 300 MS patients, patients' demographic information in summarized in Table 1. Each patient's data contain three modalities: 1) the neuroimaging data, 2) structured EHR data, and 3) clinical notes. The neuroimaging data is stored in NIFTI format. Most patients have multiple clinical visits, and during each clinic visit, a patient may have multiple laboratory tests, recorded vital signs, different prescription drugs, diagnoses, certain medical procedures, and treatments, that are recorded in the structured EHR data in separate tables. The clinical notes contain the physician's description of the patient's status at each clinic visit. Our proposed novel neural network architecture is designed to handle heterogeneous structure databases via learning representations of each modality. The prediction goal is set to a classification problem to predict if the patient has reached certain EDSS milestones at the current clinic visit. All 300 patient's EDSS scores were evaluated by physicians at the end of each clinic visit and were recorded in a table in the structured EHR. The real EDSS score was extracted from the EHR table and served as the ground truth label. The research goal is to develop a deep learning model to predict the EDSS score at the current visit using all other information and masking the true EDSS score. Figure 1 demonstrates the distributions of patients' age and EDSS. Figure 2 plots all patients EDSS historical scores along their disease course.

**Brain MRI.** We obtain a total of 360 MRI images for all 300 patients. All imaging studies were performed on a Philips 3.0T Ingenia scanner (Philips Medical Systems, Best, Netherlands). Some patients can have multiple MRIs from different clinic visits. The MRIs include five sequences: pre-contrast and post-contrast T1-weighted
sequences (T1-pre, T1-post), T2-weighted sequences, proton density-weighted sequences (PD) and fluid-attenuated inversion recovery sequences (FLAIR). All sequences were acquired with a field of view of 256 mm x 256 mm x 44 mm. For each patient, the MRI images were acquired in the axial plane. Figure 3 displays the MRI sequences of a sample patient. All MRI sequences are skull-stripped using Simple Skull Stripping (S3) [26] and the SRI24 template [27], bias-corrected using N4 Bias Field Correction to adjust the low-frequency intensity [28], and co-registered using FreeSurfer [29] to a common template (SRI24).

**Clinical Notes.** Patient's clinical notes are in free text format and contain the physician's description of the patient's health status, patient basic health information such as weight, height, BMI (body mass index), physiological status, diagnosis, medications and received treatments. We de-identified all clinical notes data by removing patients' and Physician's personal data.

**Structured EHR.** Patient's structured EHR contains laboratory tests measurements (float), vital sign observations (float), medication administrations (0/1 indicator - taken or not taken), and demographic information (age: float, race/ethnicity/gender: 0/1) in a tabular format. We construct the tables in the format of rows being observational time stamps and columns representing a number of features. The features in each table are fixed for all patients, and the number of rows for different tables and different patients is different depending on how many observational time points a patient has. For each patient's laboratory test, vital signs,s and medication table, we set the time granularity to be 4 hours. We

\begin{table}
\begin{tabular}{|l|c|c|c|c|c|} \hline  & Average \(\pm\) SD & Minimum & Maximum &.25 quantile &.75 quantile \\ \hline Age & 43.62 \(\pm\) 11.20 & 19.00 & 71.00 & 34.00 & 52.00 \\ \hline EDSS \# baseline & 1.93 \(\pm\) 1.59 & 0 & 7.50 & 1.00 & 2.50 \\ \hline EDSS \# last visit & 2.90 \(\pm\) 1.96 & 0 & 9.50 & 1.50 & 3.50 \\ \hline Number of visits & 3.39 \(\pm\) 1.60 & 1 & 13 & 2.00 & 4.00 \\ \hline Years b/w first and last visits & 5.14 \(\pm\) 4.34 & 0 & 22.66 & 2.03 & 7.01 \\ \hline Number of MRI sessions/patient & 1.20 \(\pm\) 0.96 & 0 & 4 & 0 & 2 \\ \hline \end{tabular}
\end{table}
Table 1: An overview of patient statistics in the dataset (SD: standard deviation).

Figure 1: The histograms of all patients by Age; Baseline EDSS (at initial hospital visit); EDSS at the last hospital visit; Total hospital visits; Years between the first and the last hospital visit; Number of hospital visits during which brain MRI scan was performed.
record each feature's average value if it has multiple observations during a 4-hour window. This helps to reduce table dimensions, eliminate data observational noises, and avoid creating large and sparse tables which impedes neural network training. If certain features have no values during the four-hour windows, the corresponding entry will be zero. The window size is treated as a hyper-parameter to be optimized to strike a balance between table dimension and information loss (a large averaging window smooths feature observations and blurs useful information, and a small window makes the table have a high dimension on the time-axis which and decrease the efficiency of network training). The optimal window size depends on the density of observations which can be varying for different datasets, and we found the 4-hour window is a suitable value for our dataset. We fill the entry with zero if there is no observation during that 4-hour window. Every 4-hour window is taken within a clinic encounter and can not cross two different encounters to ensure the feature values from different encounters will not be averaged together. For instance, a patient having 2 clinic encounters from 2014-05-05 1:15:00PM to 2014-05-05 6:00:00PM and 2015-09-20 9:12:00AM to 2015-09-20 1:00:00PM will have 4 rows in each table, representing the observations from 2014-05-05 12PM to 2014-05-05 4PM, 2014-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2015-05-05 12PM to 2
05 4PM to 2014-05-05 8PM, 2015-09-20 8AM to 2015-09-20 12PM, and 2015-09-20 12PM to 2015-09-20 4PM. Furthermore, we delete a row if it contains all zeros (no observations for any feature). Table 2 shows the variables we used in our dataset. All patients' demographic data is constructed as a fixed-size vector.

Figure 4 demonstrates an example patient's three clinic encounters. Note that not all data modalities were observed in each encounter.

### An encoder-decoder architecture for data fusion

We propose a multimodal neural network that takes data in various modalities (structured EHR, clinical notes, MRIs) as input and is trained to predict a patient's EDSS score. The proposed neural network adopts an encoder-decoder schema in a sequential structure with the self-attention module.

#### Encoder Network

The goal of the encoder network is to process data of different modalities and transform them into dense embeddings belonging to the same high-dimensional

\begin{table}
\begin{tabular}{c c c c c} \hline \hline  & LABORORORY TEST & \multicolumn{2}{c}{VITAL, SIGN} & MEDICATION \\ \hline Mean Composite Homoglobin & Carbon Dioxide & Albumin & Diastolic Blood Pressure & Euclidean \\ End End Clarification With & Biosynthesis & Glucose Land & Synthetic Blood Pressure & Galap
latent space. For different modalities, the neural network employs different encoder neural network architectures that are suitable for the learning task (CNN for image processing and structured EHR, Graph Neural Network or BERT for clinical notes), see Figure 5.

**Structured EHR.** The encoder network for structured EHR is composed of multiple parallel channels, and different channels embrace homogeneous network structure but with various hyper-parameters to fit the patient's structured EHR table of different sizes, see Figure 6.

The number of rows for tables can be different across different patients and across different tables of the same patient (if a certain table has no observations in a certain 4-hour window). Such an irregular-sampling format causes heterogeneity along the time axis. Traditional imputation-based methods usually define a shared regular-spaced time axis for all patients and all tables and tries to impute the missing values on the unobserved time points using various techniques, such as filling with the value of zero, average, or majority, etc., forward or backward filling, or advanced techniques such as multiple imputation [30], or the Gaussian process [3]. Imputation-based methods can be computationally expensive, time-consuming, and most importantly, increase the time dimensions, especially when the time between two clinic encounters is long.

We introduce a self-attention module into the encoder network to handle the irregularity sampling issue, by promoting the neural network to automatically pay distinctive attention (by assigning different attention weights) to different time points of patient history, aggregating them and producing a vectored representation (embedding) of each table. An attention weight is computed for each row (4-hour time window) by applying multiple layers of 1-dimensional convolutional neural networks (CNNs) on the feature dimension and outputting an attention weight for each time stamp. The attention vector can be seen as attention weights on different time stamps, and after being applied to the original input data, the network will generate an embedding of a fixed dimension that is consistent for all patients.

Figure 5: The encoder network for our proposed deep neural network.
To be specific, all channels consist of multiple stacked 1D convolution layers followed by the ReLU activation layer and dropout layers. The number of layers is set up differently for different channels according to the number of features in the input tables. For the \(i\)-th patient, the \(k\)-th data table \(\mathbf{D}_{k}^{i}\) of dimension \(t_{k}^{i}\times f_{k}\) is fed into the \(k\)-th channel, where \(t_{k}\) rows represent the time stamps of clinic visits and \(f_{k}\) columns represent variables. Note that different EHR tables (laboratory tests, vital signs, medications, etc.) have different \(f_{k}\) and different patients have different numbers of clinic visits \(t_{k}^{i}\). Each row of the table is processed through a stack of multiple 1D CNNs (see Figure 6) and is reduced to a single value (attention weight). The entire table will generate an attention weight vector \(\mathbf{\alpha}_{k}^{i}\) of size \(t_{k}^{i}\times 1\). The attention weights can be viewed as the weight factor of all \(f_{k}\) features at different time points. In the following, we omit the patient index \(i\).

We multiply the attention vector \(\mathbf{\alpha}_{k}\) with the input matrix \(\mathbf{D}_{k}\) to get the feature map \(\mathbf{e}_{k}\) for each table,

\[\mathbf{e}_{k}=\mathbf{\alpha}_{k}^{T}\cdot\mathbf{D}_{k}. \tag{1}\]

where \(\mathbf{e}_{k}\) is of size \(1\times f_{k}\). Specifically, each element in \(\mathbf{e}_{k}\) is calculated as

\[\mathbf{e}_{k}[j]=\sum_{m=1}^{t_{k}}\mathbf{\alpha}_{k}[m]\mathbf{D}_{k}[m,:],\text{ for }j=1,\ldots,f_{k}, \tag{2}\]

and \(\mathbf{e}_{k}\) is the embedding vector of the \(k\)-th table for a certain patient.

**Image Embedding.** The encoder channel for patient MRI images takes a different network structure from the structured EHR. We leverage the ResNet [31] to process the MRI images. Each MRI sequence (T1-pre, T1-post, T2, PD, FLAIR) is fed into a respective ResNet model. The output is a fixed-length latent representation of each MRI sequence as an embedding vector of a fixed preset dimension. Alternatively, other ResNet variants [32] could also serve as embedding learning network in our task. Our experiment shows the adoption of different network structure for the MRI sequences brings trivial accuracy improvements on the final prediction performance, due to the reason that 1) the ResNet model itself is powerful enough to capture the key features in the MRI images and generate diverse embeddings

Figure 6: The detailed architecture of one of the encoder channels for processing structured EHR data. The figure shows the lab test channel as an example.
for positive and negative patients; 2) the MRI data only accounts for a partial of all the input multimodal data, therefore, the effect of the ResNet variations on the final outcome will be diluted by other data modalities.

**Clinical Notes Embedding.** The encoder channel for patient clinical notes data is processed using a graph attention convolution model, which takes text as input and outputs an embedding for each document [33]. The medical word embeddings are from a pre-trained database which was trained on PubMed+MIMIC-III [34]. The graph attention model treats the entire document as a word co-occurrence network by representing words in the corpus of all patients' documents as graph nodes. In addition, we add another "document node" which represents the entire document and connects to all other nodes. The model maintains a sliding window to capture word co-occurrences, which will be represented as edges of the graph. The edge is directed and weighted, in order to represent the correct word orders in the sliding window and retain maintain meaningful semantics and word co-occurrence counts. The entire network is trained through message passing. We define \(G(V,E)\) as the graphical network, and denote a node \(v(\in V)\)'s neighbors as \(\mathcal{N}(v)\). A node \(v\) constructs a broadcasting message by aggregating (using multi-layer perceptron) its neighbor node embeddings,

\[\mathbf{m}_{v}^{t+1}=\text{ AGGREGATE }^{t+1}\left(\left\{\mathbf{h}_{w}^{t}\mid w\in \mathcal{N}(v)\right\}\right), \tag{3}\]

which can proceed in a parallel manner using matrix format,

\[\mathbf{M}^{t+1}=\text{MLP}^{t+1}\left(\mathbf{D}^{-1}\mathbf{A}\mathbf{H}^{t}\right), \tag{4}\]

where \(\mathbf{H}^{t}\in\mathbf{R}^{n\times d}\) is the \(d\)-dimensional node features of \(n\) nodes and \(\mathbf{A}\in\mathbf{R}^{n\times n}\) is the adjacency matrix, and MLP is multiple layer perceptrons neural network.

All nodes update themselves by their own embedding and all messages from their neighbors using a Gated Recurrent Unit (GRU) network,

\[\mathbf{h}_{v}^{t+1}=\text{ COMBINE }^{t+1}\left(\mathbf{h}_{v}^{t},\mathbf{m}_{v}^{t+1} \right), \tag{5}\]

again in matrix format,

\[\mathbf{H}^{t+1}=\text{GRU}\left(\mathbf{H}^{t},\mathbf{M}^{t+1}\right). \tag{6}\]

After \(T\) steps, a final self-attention read-out layer is used to aggregate all nodes embeddings and output a latent vector to represent the entire document,

\[\mathbf{Y}^{T} =\tanh\left(\hat{\mathbf{H}}^{T}\mathbf{W}_{A}^{T}\right) \tag{7}\] \[\mathbf{\beta}_{i}^{T} =\frac{\exp\left(\mathbf{Y}_{i}^{T}\cdot\mathbf{v}^{T}\right)}{\sum_{j=1} ^{n-1}\exp\left(\mathbf{Y}_{j}^{T}\cdot\mathbf{v}^{T}\right)}\] (8) \[\mathbf{u}^{T} =\sum_{i=1}^{n-1}\mathbf{\beta}_{i}^{T}\hat{\mathbf{H}}_{i}^{T} \tag{9}\]
where \(\hat{\mathbf{H}}^{T}\in\mathbf{R}^{n\times d}\) is the final node representation of all \(n-1\) nodes (remove the document node) after \(T\) time steps, and \(\mathbf{W}_{A}^{T}\) is the network parameters (a dense layer). Therefore, \(\mathbf{u}^{T}\in\mathbf{R}^{d}\) would be the final representation of the document, i.e. aggregation of all node features, which will be fed into a classification layer for document classification.

### Multi-modality Fusion

Medical data often have multiple types of information (demographics, vital, lab, diagnosis, procedure, medication, etc.) and there is intrinsic logic behind them. For example, vitals, and labs contribute to the diagnosis, and diagnosis will determine the procedure and medication. Some of this information is temporally invariant (e.g., demographics) and others are changing over time. Therefore, they need to be handled differently. Based on the causal relationships (vital, lab, MRI scan) \(\rightarrow\) (diagnosis) \(\rightarrow\) (prescription, procedure) \(\rightarrow\) (medicine administration), we build our data fusion pipeline for time-variant information through the aforementioned bidirectional GRU-based decoder. The order of the inputs (see the left part of Figure 7) are organized in a way to learn the intrinsic relationships of such information. The latent representation vectors from each encoder network channel are stacked into a regular matrix \(\mathbf{E}\) (zero-padded if not the same length), where each row represents a modality

\[\mathbf{E}=[\text{ZeroPadding}(\mathbf{e}_{1})^{T},\ldots,\text{ZeroPadding}(\mathbf{e}_ {K})^{T}]^{T}, \tag{10}\]

where \(\mathbf{E}\) is of dimension \(K\times d,d=\max(f_{1},\ldots,f_{K})\).

We integrate the time-invariant demographics at the end of the layer as a late fusion step (see the right part of Figure 7) to combine both pieces of information holistically.

### Decoder Network

We propose a decoder network structure that is composed of a stacked bidirectional GRU (Bi-GRU) network with a self-attention module taking the feature matrix \(\mathbf{E}\) as input. The self-attention serves to learn important weights on the state vectors from different data modalities. The Bi-GRU network takes \(K\) as the sequence length and \(d\) as the input size. We use \(\mathbf{C}\) to denote the stack of hidden states of all time points, which is of dimension \(K\times h,h=2\times\text{hiddensize}\) (note that factor 2 comes from the bi-direction network being used).

Each state of the bidirectional GRU network is fed into an attention module, which is 1D convolution layer of multiple output channels. The attention module outputs a vector of attention weights \(\mathbf{\gamma}\) of length \(g\) (hyper-parameter, depending on the output channel of the convolution layer), and

\[\mathbf{B}=[\mathbf{\gamma}_{1}^{T},\ldots,\mathbf{\gamma}_{K}^{T}]^{T}, \tag{11}\]

where \(\mathbf{B}\) is of dimension \(K\times g\) denoting the attention matrix. The attention matrix is multiplied with the GRU output,

\[\mathbf{O}=\mathbf{B}^{T}\cdot\mathbf{C}. \tag{12}\]
where \(\mathbf{O}\) is of dimension \(g\times h\). Note that the purpose of this attention layer is to enforce a feature reduction from the high-dimensional GRU outputs to a smaller and more informative lower-dimensional embedding not only for reducing the noise but also for increasing the efficiency of neural network training.

The output matrix \(\mathbf{O}\) is flattened, and concatenated with the patient demographic data vector \(\mathbf{d}\), and fed into a fully-connected (FC) layer for prediction,

\[o=\text{FC}(\text{Concat}(\text{Flatten}(\mathbf{O}),\mathbf{d})). \tag{13}\]

see Figure 7.

## Experiments

### Images

We introduce five channels to process the MRI sequence, where each channel employs a ResNet structure. The five channels are independent and each is trained to learn from one sequence (T1-pre, T1-post, T2, FLAIR and PD). All MRI images are bias-corrected, skull-stripped, and registered and the intensity scale is normalized [35]. If a patient performed MRI scans in more than one hospital visit, we use the last scan as it represents the patient's most recent disease status. Due to the relatively high imbalance of the positive and negative samples, we performed 10-fold re-sampling for the negative training samples during model training. To add robustness to the learning of ResNet, all input MRI images are also randomly rotated with a probability of 0.5 by a maximum of \(\pm 0.02\) degrees on all three dimensions. For each channel, a respective ResNet model is trained on the training dataset, and we select the trained model with the best performance on the validation dataset. Our goal is to learn a latent vector representation instead of performing disease classification, therefore, the training process is formulated as a metric learning task where each channel's ResNet is trained to learn an embedding for each MRI sequence of a patient. The triplet margin loss [36] operates directly on embedding distances by

Figure 7: The decoder network for our proposed deep neural network.
promoting the matching point (positive) to the reference point (anchor) and the non-matching point (negative) away from the anchor. By using a triplet margin loss, the network learns well-separated embedding vectors for positive and negative patients for downstream decoding networks to perform classification. The triplet margin loss is defined as

\[loss=\sum_{a_{i},n_{i},p_{i}\in\mathrm{batch}}\max\left(d\left(a_{i},p_{i} \right)-d\left(a_{i},n_{i}\right)+\mathrm{margin},0\right), \tag{14}\]

where \(a_{i},p_{i}\), and \(n_{i}\) are an anchor, positive and negative sample in the batch, respectively. We set the anchor point in our model as a fixed point in the embedding space, therefore, the distance from the positive samples to the anchor is minimized and the distance from the negative samples to the anchor is maximized.

The margin in the triplet margin loss is chosen to be 1.5. The learning rate is set to be \(10^{-5}\) and the batch size is 10. The ResNet in each encoder channel is trained for 500 epochs. Early stopping criteria of not-improving for consecutive 50 epochs on the validation dataset is adopted.

We leverage the gradient-weighted class activation mapping (Grad-CAM) [37] model to locate and visualize the important regions the ResNet neural network is learning for predicting the target. The Grad-CAM uses flowing gradients of the prediction target into the last convolutional layer of the ResNet to produce a heatmap of the regions according to their contributions to the prediction, see Figure 8,

Figure 8: Attention maps for MRI sequences of a sample patient.
### Clinical Notes

We preprocess patients' clinical notes by identifying and then removing all sensitive patient health information that is irrelevant to our prediction task, including the patient and physician's name, address, phone number, and email address. Similar to the MRI image data, we formulate the embedding generation problem from clinical notes as a metric learning problem, where the message-passing graph neural network is trained to learn meaningful embeddings and their distances between positive and negative samples. Hence, the same loss function (14) is used for this encoder channel.

We set the size of the window to be 10 (covering 10 consecutive words) and the message passing layer to be 2. The hidden side of the GRU network is 64. We trained the graph network with 500 epochs with a batch size of 128, the learning rate of \(10^{-3}\), and early stopping criteria of 50 epochs (no improvements on the validation dataset). We choose the best-performing model on the validation dataset and run it on the test dataset to get the model's final performance.

### Structured EHR

The patient's structured EHR consists of tables of 4 categories, the laboratory tests table, the vital signs table, the medications table, and the demographics table. The first 3 categories are in the format _number of timestamps \(\times\) number of features_ containing the laboratory test results (float), vital sign measurements (float), and medications (0/1 indicators), respectively. Table 2 shows a pre-selected subset of all the variables from the above 3 categories to be used in our model, based on their observation frequency. The demographic table contains race (0/1, one-hot encoded), ethnicity (0/1, one-hot encoded), sex (0/1, male/female), and age (float, min-max normalized). The encoder network consistents 3 channels for each of the first 3 categories, and the network parameters are described in Table 3.

A patient's three structured EHR's embeddings produced by the encoder network will be concatenated with the five MRI image embeddings produced by the ResNet, and together with the clinical note embedding to be fed into the decoder network. In the situation of a patient (a small amount) without MRI or clinical notes, the corresponding embedding will be set to an all-zero vector. In the decoder network, the bidirectional GRU network is set to have 4 layers and hidden size of 512.

### Results

We used 5-fold cross-validation by randomly split the 300 patients into five folds and iteratively using each fold as the hold-out test set (20%) and the remaining as the training set (80%). The model's performance of predicting EDSS \(>\) 4.0 using different data modalities and their combinations are presented in Table 4. The prediction goal is whether the patient's EDSS \(>\) 4 of the current clinic visit, using longitudinal data of the current clinic visit and all previous visits. Multimodal data inputs in general perform much better than single modal input, and the top-3

\begin{table}
\begin{tabular}{l c c c c c c c} \hline \hline  & Candid & Dupont & Candid & RLO \(>\) Dupont & Comid & RLO \(>\) Dupont & Pougles \\ \hline Observed 1 (left) & 1.0 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & Age \\ Observed 2 (Vital Sign Observ.) & 1.1 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & Age \\ Observed 3 (Medicine) & 1.1 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & Age \\ \hline \hline \end{tabular}
\end{table}
Table 3: Encoder network parameters (I: input channel size, O: output channel size, K: kernel size, S: stride size, P: padding size, R: (dropout) rate).
AUROC performances are when all data (0.9301), EHR & Notes (0.9196), and MRIs & Notes (0.9201) are being used. The degradation in performance by deleting MRI or EHR information from the input data is very limited. However, if clinical notes were to be deleted, the performance drops to 0.8499. Table 5 shows the model's performance for predicting other EDSS milestones (EDSS\(>\)6 and EDSS\(>\)7) using all data modalities.

The encoder channels (laboratory tests, vital signs, medications) for the structured EHR, as a self-attention network by itself, could also generate the feature importance. Notice that the feature importance can be generated both on an individual-level and global level. The latter is evaluated as the average of all feature importance of all individuals. Figure 9 shows the global importance of all laboratory features evaluated on all patients in the test set, where a larger value corresponds to the higher importance of a feature. The top 3 important features for all patients, as can be seen from the figure, are "Absolute Neutrophils", "Absolute Lymphocytes" and "Absolute Monocytes". Figure 10 shows the global-feature importance on all vital signs and medications, some medicine such as the "Baclofen 10 MG Oral Tablet", "Gabapentin 300 MG Oral Capsule", "predniSONE 50 MG Oral Tablet" were commonly used to treat MS symptoms which were identified and assigned with high importance by our algorithm. For the attention weights on vital signs, the "Temperature", "Respiration", "Pulse Quality" and "Respiration Quality" are reasonably assigned with the least importance for the prediction of MS.
brain MRI images and the clinical notes, and structured EHR contributes the least to this prediction goal.

The results shows that despite the many publications, conventional MRI contains relatively less information about MS severity compared to other data modalities. However, T2 and Flair MRI performed relatively better than other MRI sequences. Clinical notes were well-documented to be used for the prediction of EDSS, which has been re-verified in our experiment as the relatively not good performance of using MRIs, or EHR, or MRIs & EHR were all improved when clinical notes were added to the input. A re-examination of the data reveals a reasonable explanation that the clinical notes contain rich patient general disease information including patient status, medical procedures, and treatment information, which implicitly and partially embeds information from the EHR data and MRI images.

Future research directions and limitations

The study focuses on predicting patients' MS severity as current clinic visits by using current and historical medical information, with the goal to develop an AI

Figure 9: The attention weights for laboratory tests.
based patient disease status evaluation tool to replace the human expert. A more interesting research question is to predict a patient's MS disease progression in the future. This necessarily needs to consider the EDSS change "rate" by considering the disease duration. For example, having an EDSS of 4.0 at age of 65 and a disease duration of 40 years would mean a relatively benign disease but having an EDSS score of 4.0 only after 5 years of MS diagnosis is considered as "aggressive" MS. Moreover, as one of the reviewers that has helped improve this paper pointed out, the severity of MS can be seen as a relative concept instead of an absolute one. The severity of MS should be studied based on an understanding of the "natural" disease progression, and it varies in terms of many factors (eg. sex, disease duration, lesion load, atrophy, etc.) Limited by the data size and commonly agreed on criteria to distinguish the "aggressive" cases from the rest, we focus on developing a tool to predict EDSS milestones at the moment and leave the decision of MS severity to MS specialists by jointly considering all the above factors. In addition, this problem itself is quite an interesting research problem and could potentially be studied using survival methods, the results will have a high impact on the prevention of rapid disease progression through early intervention.

The second is the limitation of the imaging data. While random rotation of MRI scans (a data augmentation technique used to train ResNet on the MRI sequences) helps generalizability, the use of only one scanner for all datasets makes it difficult to infer if the model would work in the same way when introduced to new images from a different scanner. Therefore, our work serves as a proof-of-concept regarding this question. Ideally, more data (especially data from external sources) needs to be collaboratively collected to verify the inclusion of MRI potentially has a positive impact on a multi-modal model.

Thirdly, the study was conducted on a cohort of 300 MS patients from a local academic medical center. An important future research direction is to evaluate the generalizability of the proposed model to other institutions. The result replicability should be checked from two perspectives, the first is the prediction accuracy with or without model retraining, and the second is if the ranking of importance for different data modalities is the same in general, for example, MRI images and clinical notes contains more signals compared to the structured EHR. If the results in this study are verified, it may serve as a cost-effective study recommending which electronic health information should be collected to reach maximum prediction accuracy.

Figure 10: The attention weights for (a) medications and (b) vital signs.
## Abbreviations

MS: Multiple sclerosis; MSSS: Multiple sclerosis severity scores; P-MSSS: Patient-derived MSSS; EDSS: Expanded disability status scale; EHR: Electronic health records; AUROC: Area Under the receiver operating characteristic curve; AUPRC: Area under the precision-recall curve MRI: Magnetic resonance imaging; SD: standard deviation; PD: Proton density; GRU: Gated Recurrent unit; Grad-CAM: Gradient-weighted class activation mapping; BCE: binary cross entropy.

## Supplementary Information

### Acknowledgements

We thank the reviewers for proposing a critical perspective of viewing this research problem which greatly helped improve the quality of this manuscript.

### Authors' contributions

KZ, S5, and XJ conceived the original idea and designed the model. KZ implemented the model and conducted the experiments. EB and JL contributed to the conception, data acquisition, and interpretation. KZ formed the manuscript. SE, EB, JL, and XJ critically revised the manuscript. EB and JL were in charge of overall direction and planning and helped supervise the project.

### Funding

XJ is CPRIT Scholar in Cancer Research (RR180012), and he was supported in part by Christopher Saroffin Family Professorship, UT Stars award, UTHealth startup, the National Institute of Health (NIH) under award number R01AG066749 and U01TR002062,

### Availability of data and materials

The data that support the findings of this study are available on request from the corresponding author SS. The data are not publicly available due to their containing information that could compromise the privacy of research participants. Code is publicly available on Github: [https://github.com/anotherkaihang/MS](https://github.com/anotherkaihang/MS).

### Declarations

**Ethical approval and consent to participate**

The study protocol was approved before the initiation of this study by the Committee for the Protection of Human Subjects of the University of Texas Health Science Center at Houston under IRB: HSC-MS-02-090. All recruited patients provided written informed consent upon enrollment. All methods were performed in accordance with the Declarations of Helsinki.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that there is no conflict of interest.

### Author details

1School of Biomedical Informatics, University of Texas Health Sciences Center at Houston, Houston, TX, United States. 2Department of Neurology, University of Texas Health Sciences Center, McGowern Medical School, Houston, TX, United States. 3Division of General Internal Medicine, Department of Internal Medicine, University of Texas Health Sciences Center, McGowern Medical School, Houston, TX, United States. 4Department of Applied Data Science, San Jose State University, San Jose, CA, United States.
