Phenotypic Diversity in Multiple Sclerosis Can Be Represented by Four Additive Symptom Modules

1. Introduction

Multiple sclerosis (MS), unlike many neurological diseases, does not have a single recognizable phenotypic core. Alzheimer’s disease is typically anchored by memory impairment, Parkinson’s disease by bradykinesia, rigidity, and resting tremor, and Charcot–Marie–Tooth disease by sensory loss, hyporeflexia, weakness, and muscle atrophy. In MS, by contrast, some patients present with sensory symptoms [1], weakness [2], optic neuritis and visual loss [3,4], or ataxia and incoordination [5,6,7]. Other presentations of MS include impaired cognition [8], psychological disturbances [9], fatigue [10,11], spasticity [12], pain [13,14], gait disorders [15], sphincter disturbances [16], and paresthesias [17]. No single symptom complex defines all cases at onset or during the course of the disease. Traditional classifications by disease course (relapsing–remitting, secondary progressive, primary progressive) capture the temporal pattern of MS but do not fully describe the cross-sectional distribution of symptoms across functional systems.

Disease progression in MS is also more complex than simple worsening along a single clinical axis. In diseases such as Alzheimer’s disease and Parkinson’s disease, progression often consists of gradual worsening of the core symptoms that define the disorder at onset. In MS, initial symptoms may worsen over time, but new episodes may also add deficits in domains that were not previously affected. Thus, clinical worsening may reflect both increasing severity within an existing symptom domain and expansion into new domains.

As MS evolves, some patients remain relatively restricted to one or two domains of impairment, whereas others accumulate deficits across multiple domains.[18,19] Kister et al. [20] examined symptom expansion in more than 25,000 MS patients across 11 symptom domains, including mobility, hand function, cognition, fatigue, bowel-bladder function, sensory symptoms, spasticity, pain, depression, and ataxia-tremor. They found a steady expansion in symptom extent with increasing disease duration. This pattern suggests that MS progression reflects the additive recruitment of multiple symptom domains rather than progressive worsening of a single invariant symptom core. Although relatively pure phenotype profiles may be observed, most patients can be viewed as reflecting admixtures of phenotype features that vary in both severity and domain.

Kurtzke recognized this modular structure clinically and designed the Functional System Scores (FSS) to rate relatively independent neurological systems, including pyramidal, cerebellar, brainstem, sensory, visual, bowel and bladder, and cerebral functions. Subsequent analyses have shown that these systems are only weakly to moderately correlated, supporting their interpretation as partially independent axes of impairment rather than as indicators of a single unidimensional disability construct. Although Hobart and colleagues [21] argued that the FSS should not be summed to yield a psychometric total disability score, the pattern of low–moderate intercorrelations is consistent with the view that MS phenotype is composed of separable symptom modules. Such a multi-axis representation is naturally amenable to decomposition by non-negative matrix factorization.

Non-negative matrix factorization [22,23] (NMF) is a linear multivariate method that decomposes a non-negative data matrix X into the product of two lower-rank non-negative matrices W and H, such that $X\approx WH$. The rows of W represent patients as mixtures of latent modules, and the rows of H specify the contribution of each observed feature to those modules. Because all entries are constrained to be non-negative, NMF represents each patient as an additive combination of modules rather than as a mixture of opposing positive and negative loadings. In practical terms, NMF reduces a larger set of observed features to a smaller number of additive modules that summarize recurrent patterns in the data.

In exploratory analysis of high-dimensional patient phenotype data, an important goal is to obtain a low-dimensional representation that supports both visualization and interpretation. Traditional unsupervised workflows often cluster patients in a high-dimensional feature space and then visualize them with nonlinear embedding methods such as t-SNE or UMAP. Although these approaches can reveal regional structure, the resulting axes are not directly interpretable and the clinical meaning of cluster membership may be difficult to recover. In contrast, NMF directly reduces correlated symptom features to additive modules and represents each patient as a weighted mixture of those modules. This additive constraint distinguishes NMF from methods such as principal component analysis (PCA) or independent component analysis (ICA), which allow negative loadings and therefore represent data as combinations of opposing components.

In neurological disease, this additive structure is not merely a mathematical convenience but may reflect clinical reality. In MS, deficits accumulate over time as lesions affect distinct motor, sensory, autonomic, and cognitive systems. The observed phenotype can therefore be viewed as a superposition of deficits across multiple functional domains rather than as movement along a single severity axis. Under this interpretation, NMF provides a natural framework for modeling MS phenotype as a weighted combination of latent symptom modules.

We therefore asked whether the known heterogeneity of MS phenotype, as documented in routine clinical notes, could be summarized by a small number of additive symptom modules and whether individual patients could be represented as relatively pure or admixed combinations of those modules. The purpose of this approach was not simply to reduce 17 phenotype features to four modules, but to create a quantitative phenotype space in which each patient can be described by both a predominant clinical pattern and an effective number of contributing modules.

2. Methods

2.1. Neurology Notes

A total of 12,661 de-identified neurology clinical notes were obtained from a REDCap database at the University of Illinois Hospital, the primary teaching hospital of UI Health, for the period from January 14, 2016, through September 2, 2022. Notes were eligible for inclusion if they were longer than 600 words, associated with a diagnosis of multiple sclerosis (ICD-10-CM G35), and classified as Progress Note encounters. Patient-directed Visit Summaries were excluded. After deduplication, 4,617 notes remained for analysis. Use of de-identified clinical documentation for research was approved by the Institutional Review Board of the University of Illinois (Protocol No. 2017-0520Z, amendment approved August 10, 2022).

Each note represented a narrative summary of a patient encounter and served as the unit of analysis. Prior to annotation, notes were converted to JavaScript Object Notation (JSON) format, with each note represented as a single JSON object. The primary data field in each object was note_text, which contained the full note in plain-text format. Associated metadata included note_id, patient_id (de-identified), age, sex, race, note_length (characters), note_count (number of notes available for that patient), note_date, note_type (e.g., Progress Note), clinic (e.g., Neurology), and elapsed_days (days from the first note to the current note) (see Appendix Table A1).

2.2. Categorization of Neurology Notes by Phenotype

Each neurology note was categorized with respect to 17 non-mutually exclusive neurological phenotype features, including weakness, sensory symptoms, pain, ataxia, cognitive impairment, bladder symptoms, bowel symptoms, fatigue, gait impairment, and visual symptoms. Phenotyping was performed using the OpenAI API and GPT-5.2, following previously described and validated note-level methods [24]. Briefly, the model was prompted to classify each complete clinical note for the presence or absence of each phenotype feature, yielding a structured binary note-level feature vector. In the validation study, 100 de-identified MS neurology progress notes were independently annotated by two human annotators, discordant labels were adjudicated, and GPT-5.2 was evaluated against the adjudicated reference set. GPT-5.2 achieved macro-precision of 0.734, macro-recall of 0.921, macro-F1 of 0.801, and macro-averaged Matthews correlation coefficient of 0.777, approaching the performance of human annotators. The complete label definitions and prompt structure have been reported previously [24]; the 17 phenotype features and definitions used in the present analysis are summarized in Appendix Table A1.

Because most patients contributed multiple notes over time, note-level phenotype assignments were aggregated at the patient level. For each patient, the full set of notes was summarized as a 17-dimensional binary phenotype burden vector, where 1 indicated that a phenotype feature was present in at least one note and 0 indicated that the feature was absent from all notes. Repeated mentions of the same phenotype within a single note or across multiple notes did not increase the score. After aggregation, patients with no phenotype features were excluded, leaving 578 patients for non-negative matrix factorization. The final phenotype burden vector represented cumulative phenotype burden across all available notes, with each dimension constrained to binary values.

2.3. Non-Negative Matrix Factorization

Non-negative matrix factorization (NMF) was performed using the NMF implementation in sklearn.decomposition, with max_iter = 5000, random_state = 0, and init = nndsvd. We examined 2-, 3-, 4-, and 5-module solutions by varying n_components.

For each candidate rank, model fit was summarized using relative reconstruction error, computed as the Frobenius norm of the residual matrix divided by the Frobenius norm of the original data matrix:

$${\displaystyle \frac{{\parallel X-WH\parallel }_F}{{\parallel X\parallel }_F}}$$

(1)

Here, X denotes the original patient-by-feature phenotype matrix, $WH$ denotes its NMF reconstruction, and ${\parallel ·\parallel }_F$ denotes the Frobenius norm, equivalent to the square root of the sum of squared matrix entries. Lower relative reconstruction error indicates better reconstruction of the original phenotype matrix.

Selection of the number of NMF components is a recognized challenge in unsupervised EHR phenotyping; prior work has emphasized that model order should be chosen by balancing reconstruction error, stability, interpretability, and clinical validity rather than by fit alone [25].

Using the highest-loading feature as a provisional label, heat maps were constructed for the preferred 4-module solution (Figure 2) and for the competing 2-, 3-, and 5-module solutions (Figure A1–Figure A3). Although the 4-module solution was intermediate with respect to reconstruction error (Table 1), it showed the greatest clinical interpretability and was selected as the preferred solution. The four modules were named sensory-visual-pain, ataxic-spastic-falls, cognitive-psychologic-fatigue, and autonomic-bladder-bowel on the basis of their three highest-loading phenotype features (Table 2).

To assess stability of the selected four-module solution, we performed 100 repeated 90% subsampling runs. In each run, a random 90% sample of patients was selected without replacement, the four-module NMF model was refit, and the resulting module-loading vectors were matched to the corresponding full-cohort modules. Because NMF components are unordered, modules were matched using one-to-one maximum-correlation assignment based on Pearson correlation of the 17-dimensional feature-loading vectors. Stability was summarized as the median, interquartile range, minimum, and maximum matched correlation for each module (Appendix Table A3).

2.4. Computed Features

Using the patient-level weights from the four-module NMF model, the contribution of each module was normalized to sum to 100% for each patient. The dominant module was defined as the module with the highest normalized loading. Patients were classified as highly dominant if one module accounted for more than 55% of the total loading and as near-pure if one module accounted for more than 70% of the total loading. Patients with 100% of their normalized loading in a single module were classified as single-module, or pure cases.

Entropy was calculated from the row-normalized NMF weights, with each patient’s four module weights scaled to sum to 100%. The effective number of modules was then computed as $e^H$, where H is Shannon entropy.

2.5. Plotting and Visualization

Heat maps showing feature loadings for the 2-, 3-, 4-, and 5-module NMF solutions were created using the seaborn library. Because the four-module solution was selected as the preferred model, patient-level module weights were normalized to sum to 100%, allowing each patient to be represented within a tetrahedral, or 3-simplex, phenotype space. In this representation, each vertex corresponds to a pure module profile, edges correspond to two-module mixtures, faces correspond to three-module mixtures, and interior points correspond to admixtures of all four modules. The tetrahedron plot was created with Matplotlib using mpl_toolkits.mplot3d.

To provide complementary views of individual triangular faces of the tetrahedron, simplex plots were created using the python-ternary library in combination with Matplotlib. Histograms and box plots of entropy-derived effective module counts were created with Matplotlib.

3. Results

We used a large language model to extract neurological phenotype features from 4,617 neurology notes representing 578 patients with multiple sclerosis. Phenotypes were identified and categorized at the whole-note level [24]. Features were aggregated across all available notes for each patient, yielding a single binarized 17-dimensional phenotype vector per patient. The most common features were abnormal gait, pain, weakness, and abnormal sensation (Figure 1). Patients with no phenotype features were excluded. Among the 578 patients retained for analysis, the mean number of phenotype features per patient was 8.09 ± 3.35 (minimum = 1, maximum = 16, median = 8).

We applied non-negative matrix factorization to the patient phenotype vectors and selected the 4-module solution (Figure 2) as the preferred model based on reconstruction error, clinical interpretability, parsimony, and stability (Table 1; Table 2; Appendix Table A3). The four modules identified by non-negative matrix factorization—sensory-visual-pain, ataxic-spastic-falls, cognitive-psychologic-fatigue, and autonomic-bladder-bowel—were clinically interpretable and corresponded to recognizable patterns of symptom clustering in MS. The highest-weighted features within each module were neurologically coherent and aligned with plausible patterns of dysfunction across sensory and visual pathways, pyramidal-cerebellar motor systems, cognitive-affective networks, and autonomic or spinal cord pathways.

Figure 3 provides a two-dimensional summary of the four-module compositional space by separating module identity from module admixture. The x-axis identifies the dominant phenotype module, whereas the y-axis quantifies phenotypic admixture as the entropy-derived effective number of modules. Thus, each patient is summarized by both a predominant clinical direction and a quantitative measure of how broadly the phenotype is distributed across modules. Patients were classified as highly dominant if one module accounted for more than 55% of the total loading and as near-pure if one module accounted for more than 70% of the total loading. Patients with 100% of their normalized loading in a single module were classified as pure single-module, or pure-vertex, cases. The distribution of dominant, highly dominant, near-pure, and pure module assignments is summarized in Table 3.

The 4-module solution also allowed each patient to be represented as a normalized four-dimensional compositional vector, with the four module weights summing to 100%. Because all module weights were non-negative and row-normalized to sum to 100%, each patient could be visualized within a regular tetrahedral space, or 3-simplex. In this representation, each vertex corresponds to a pure module profile, edges correspond to two-module mixtures, faces correspond to three-module mixtures, and interior points correspond to admixtures of all four modules. Patients mapped near a vertex in the tetrahedral projection can therefore be interpreted as near-pure examples of a single module, whereas patients located in the interior represent more broadly admixed phenotypes. Thus, the tetrahedron should be interpreted as a convenient visualization of normalized four-module composition rather than as independent evidence for discrete patient clusters (Figure 4).

The tetrahedron has four faces, each bounded by three vertices representing pure modules. Two of the four faces are shown in Figure A4 and Figure A5.

We also examined whether the ideal tetrahedral vertices were occupied by empirical patient cases. Only seven patients occupied pure vertices, corresponding to two cognitive-psychologic-fatigue and five sensory-visual-pain profiles (Table 3). The remaining two ideal vertices, corresponding to autonomic-bladder-bowel and ataxic-spastic-falls profiles, were not occupied by pure patient cases. These findings indicate that the four-module NMF solution defines an ideal compositional phenotype space, but the observed cohort did not fully occupy all pure module vertices.

As an additional measure of phenotypic admixture, we calculated the Shannon entropy of the four normalized module weights for each patient. When 100% of the loading is concentrated in a single module, entropy is 0; when the loadings are distributed equally across all four modules, entropy reaches its maximum value of $ln\left(4\right)=1.386$. Entropy values were converted to the effective number of modules using the formula $e^H$, where H is Shannon entropy. The effective number of modules therefore ranges from 1 to 4, with values near 1 indicating that a patient’s phenotype is concentrated within a single module and values near 4 indicating maximal admixture across all four modules (Figure 3 and Figure 5). Importantly, a patient may have multiple phenotypic features within the same module group, such as sensory symptoms, pain, and visual symptoms, while still having an effective module count near 1 if those features map predominantly to a single latent module.

Figure 1. Patient counts by phenotype feature. Most common features were gait disorders, pain, and sensory loss or paresthesias. Counts are across all available notes but each feature is counted only once per patient.

Figure 1. Patient counts by phenotype feature. Most common features were gait disorders, pain, and sensory loss or paresthesias. Counts are across all available notes but each feature is counted only once per patient.

Figure 2. Heatmap of feature loadings for the preferred 4-module NMF solution with sensory-visual-pain, ataxic-spastic-falls, cognitive-psychologic-fatigue, and autonomic-bladder-bowel as primary phenotypes of each module. Note that weakness and gait load significantly on Modules M1 and M2.

Figure 2. Heatmap of feature loadings for the preferred 4-module NMF solution with sensory-visual-pain, ataxic-spastic-falls, cognitive-psychologic-fatigue, and autonomic-bladder-bowel as primary phenotypes of each module. Note that weakness and gait load significantly on Modules M1 and M2.

Figure 3. Effective number of phenotype modules by dominant NMF module. Each patient was assigned to the module with the largest normalized NMF weight. The y-axis shows the entropy-derived effective number of modules, which quantifies phenotypic admixture. Values near 1 indicate relatively pure, module-dominant profiles, whereas values near 4 indicate broad admixture across all four modules. The figure shows that most patients retain contributions from multiple modules despite having an identifiable dominant module. Boxes show the median and interquartile range; whiskers extend to 1.5 times the interquartile range. Individual points represent patients.

Figure 3. Effective number of phenotype modules by dominant NMF module. Each patient was assigned to the module with the largest normalized NMF weight. The y-axis shows the entropy-derived effective number of modules, which quantifies phenotypic admixture. Values near 1 indicate relatively pure, module-dominant profiles, whereas values near 4 indicate broad admixture across all four modules. The figure shows that most patients retain contributions from multiple modules despite having an identifiable dominant module. Boxes show the median and interquartile range; whiskers extend to 1.5 times the interquartile range. Individual points represent patients.

Figure 4. Module weightings for each patient plotted in four-dimensional space, where each vertex represents a relatively pure presentation of one module phenotype. Panel A shows patients colored by dominant module when the dominant module contributed >0.55 of the total phenotype loading. Panel B shows patients colored by dominant module when the dominant module contributed >0.70 of the total phenotype loading. Patients shown in gray have admixed module phenotypes.

Figure 4. Module weightings for each patient plotted in four-dimensional space, where each vertex represents a relatively pure presentation of one module phenotype. Panel A shows patients colored by dominant module when the dominant module contributed >0.55 of the total phenotype loading. Panel B shows patients colored by dominant module when the dominant module contributed >0.70 of the total phenotype loading. Patients shown in gray have admixed module phenotypes.

Figure 5. Histogram of the effective number of phenotype modules. The effective number of modules quantifies the degree of admixture among the four phenotype modules for each patient. It is derived from Shannon entropy and ranges from 1.0, indicating a highly pure phenotype assignment, to 4.0, indicating maximal admixture. Most patients had an effective number of modules greater than 2.5 and were considered moderately to highly admixed.

Figure 5. Histogram of the effective number of phenotype modules. The effective number of modules quantifies the degree of admixture among the four phenotype modules for each patient. It is derived from Shannon entropy and ranges from 1.0, indicating a highly pure phenotype assignment, to 4.0, indicating maximal admixture. Most patients had an effective number of modules greater than 2.5 and were considered moderately to highly admixed.

4. Discussion

In this study, we modeled the phenotype of multiple sclerosis as a non-negative mixture of additive symptom modules derived from empirical clinical-note data. The motivation for this approach is grounded in the long-recognized polysymptomatic nature of MS [6]. MS is defined as an autoimmune demyelinating disease characterized by lesions disseminated in time and space [26]. Kurtzke’s Expanded Disability Status Scale (EDSS) and Functional Systems Scores (FSS) [27,28,29] codified this clinical heterogeneity by distributing disability across pyramidal, cerebellar, brainstem, sensory, bowel and bladder, visual, and cerebral systems [29]. Although the FSS provides a clinically intuitive multi-axis framework, it assigned signs and symptoms to functional systems using clinician-defined categories. In contrast, the modules identified here were derived from data-driven co-occurrence patterns among note-derived phenotype features.

This distinction is important because MS phenotype is not stereotyped. Individual patients accumulate different combinations of sensory, motor, cognitive, and autonomic impairments over time, reflecting demyelinating lesions distributed across different CNS locations and disease episodes.

The empirical feature distributions in this cohort supported this view. Most patients had multiple phenotype features documented across their available notes, and the most frequent features included gait impairment, pain, weakness, and sensory symptoms (Figure 1). Thus, the input data did not consist of isolated single-symptom presentations, but of broad patient-level phenotype profiles accumulated over time. This polysymptomatic structure provided the rationale for modeling MS phenotype as a mixture of latent modules rather than as a set of mutually exclusive clinical categories.

The observed phenotype can be viewed as the additive superposition of system-specific impairments rather than as progression along a single severity axis. NMF provides a natural framework for representing each patient as a non-negative mixture of latent symptom modules. We examined 2-, 3-, 4-, and 5-module solutions and selected the 4-module solution because it provided the clearest clinical interpretation while maintaining reasonable reconstruction error. This choice should be viewed as a balance between model fit and interpretability rather than as evidence for a uniquely optimal number of biological subtypes.

Our four-module solution identified clinically coherent modules corresponding to sensory-visual-pain, ataxic-spastic-falls, cognitive-psychologic-fatigue, and autonomic-bladder-bowel patterns (Table 2). These empirically derived modules broadly resembled familiar neurological functional domains, including sensory/visual pathways, pyramidal-cerebellar motor systems, distributed cognitive-affective networks, and autonomic/spinal cord pathways with supraspinal modulation.

The most practical summary of this representation is two-dimensional. Each patient can be described by the identity of the dominant module and by the effective number of modules contributing to the phenotype (Figure 3). The dominant module identifies the principal clinical direction of the phenotype, whereas the effective number of modules quantifies the degree of phenotypic admixture. This distinction is important because two patients may share the same dominant module but differ substantially in whether their phenotype is relatively pure or broadly distributed across multiple modules.

When the four non-negative module weights were normalized to sum to 100%, patients could be visualized within a tetrahedral, or 3-simplex, space. This visualization provides an intuitive geometry for module composition: vertices represent pure module profiles, edges represent two-module mixtures, faces represent three-module mixtures, and interior points represent admixtures of all four modules. The tetrahedron should be interpreted as a convenient visualization of normalized four-module composition rather than as independent evidence for four discrete patient clusters. In this study, most patients occupied admixed regions of the interior phenotype space, although a subset showed relatively pure or strongly dominant patterns near the vertices. These purer phenotypes were most often sensory-visual-pain dominant, whereas relatively pure ataxic-spastic-falls, cognitive-psychologic-fatigue, and autonomic-bladder-bowel phenotypes were less common. Only seven patients occupied true vertices of the tetrahedron, corresponding to pure sensory-visual-pain and pure cognitive-psychologic-fatigue profiles (Table 3).

The presence of relatively pure vertex-adjacent phenotypes invites comparison with archetypal analysis, a method designed to identify extreme representative profiles on the convex hull of the data. Archetypal approaches have been used to define clinically interpretable extreme disease states from longitudinal clinical data. For example, Trescato et al. [30] derived archetypal phenotypes in amyotrophic lateral sclerosis and used them to model disease-progression trajectories, illustrating the value of low-dimensional phenotype representations for studying clinical heterogeneity. The present analysis differs in an important respect: whereas archetypal analysis identifies extreme representative disease states, non-negative matrix factorization identifies additive latent phenotype modules that may coexist within the same patient. Thus, the vertices of the tetrahedron should not be interpreted as archetypes discovered by archetypal analysis, but rather as limiting cases of a normalized four-module composition.

The predominance of admixed phenotypes is consistent with clinical experience: as MS evolves, patients often accumulate deficits across multiple functional systems rather than remaining confined to a single domain [20]. At the same time, the existence of relatively module-pure patients suggests that some individuals may have phenotypes dominated by particular functional systems. Whether these purer patterns reflect biologically meaningful subtypes, stochastic lesion distribution, differences in disease stage, or differences in documentation remains uncertain. If module-dominant groups show distinct MRI lesion topography, regional atrophy, immunologic signatures, fluid biomarkers, genetic risk profiles, or treatment responses, this would support the biological relevance of the module structure. Conversely, if biological markers are similar across module-dominant groups, the modules may primarily reflect the geometry of CNS functional organization under a common pathogenic process.

The principal practical value of this approach is that it converts heterogeneous clinical-note phenotypes into quantitative patient-level module scores. These scores have three components: 1) the identity of the dominant symptom module, 2) the magnitude of the dominant module, and 3) the effective number of modules. These scores may serve as predictors, outcomes, or stratification variables in future models of relapse risk, disability progression, treatment response, MRI lesion distribution, regional atrophy, or biomarker profiles.

Conventional one-dimensional severity scales do not explicitly capture both module dominance and degree of admixture. Even multi-axis clinical instruments such as the Kurtzke Functional Systems Scores provide domain-specific scores but do not directly yield a normalized measure of phenotypic admixture such as the effective number of modules. In this sense, the four-module representation has a regularization-like effect: it trades some feature-level granularity for a more stable and interpretable description of recurring phenotypic patterns.

These findings do not establish distinct MS disease entities. Rather, they provide a quantitative framework for representing MS phenotypic diversity and for generating testable hypotheses about clinically or biologically meaningful subgroups. Future studies should determine whether these NMF-derived module scores are reproducible across cohorts and whether they predict independent measures of disease activity, progression, treatment response, or biological mechanism.

4.1. Relation to Prior MS Phenotype Studies

Prior MS phenotyping studies generally begin with a rectangular patient-by-feature matrix and follow one of two complementary strategies (Table 4). One strategy uses phenotype features to group patients into a small number of clinically interpretable classes or clusters, making the patient the primary object of analysis. The other examines the structure of the phenotype features themselves, reducing correlated features into a smaller set of latent dimensions, factors, or communities. Thus, while both approaches analyze the same m × n patient-feature matrix, they differ in whether the focus is on categorizing patients or identifying underlying feature structures.

Shahrbanian et al. [31] used hierarchical clustering to examine relationships among nine MS symptom variables and identified three broad variable clusters: cognitive-emotional symptoms, pain-fatigue-sleep symptoms, and spasticity-balance symptoms. Although clinically intuitive, this analysis produced a relatively coarse symptom grouping rather than a patient-level compositional phenotype model. Gulick [32] used factor analysis to reduce 22 phenotypic features to 5 factors: skeletal-motor (weakness, ataxia, spasticity), elimination (bowel and bladder), emotional (depression and anxiety), sensory, and brainstem/visual.

Ajdacic-Gross et al. [33] examined 20 symptom phenotypes in 1942 multiple sclerosis subjects and derived six distinct classes through latent class analysis, which included multiple symptoms (14.1%), gait-balance (13.5%), fatigue-weakness (11.7%), gait-paralysis (23.9%), vision (21.3%), and paresthesia (15.4%). De Nadai et al. [34] applied latent profile analysis to 11 Multiple Sclerosis Patient-Reported Outcome (MS-PRO) scales, which captured patient-reported impairment across mobility, hand function, vision, fatigue, cognition, bladder/bowel function, sensory symptoms, spasticity, pain, depression, and tremor/coordination domains. They derived 9 distinct subtypes including normal functioning, fatigue, sensory, somatic, somatic plus cognitive, severe disability, poor mobility, moderate disability, and physical symptoms. Howlett-Prieto et al. [35] applied the Louvain algorithm to detect network communities in a 113 patient by 17 symptom feature array and found five unipartite communities (pain, fatigue, cognitive, sensory, and gait-weakness-hypertonia).

Taken together, these studies support the view that MS phenotype is multidimensional and cannot be adequately represented by a single severity axis. However, they differ from the present study in three important respects. First, several prior studies used patient-reported outcomes or symptom-onset questionnaires, whereas the present analysis used phenotype features extracted from longitudinal clinical notes. Second, clustering, latent class analysis, and latent profile analysis assign patients to discrete groups, whereas NMF represents each patient as a mixture of additive symptom modules. Third, prior cluster or class solutions generally produced mutually exclusive patient categories, while the present approach explicitly quantifies admixture through normalized module weights, module dominance, entropy, and the effective number of modules. Thus, the present findings should be viewed not as a replacement for prior MS phenotyping studies, but as a complementary representation of MS phenotype as a continuous, modular, and compositional space.

4.2. Limitations and Future Directions

This study has several limitations that also define important directions for future work. First, phenotyping was performed by a large language model rather than by manual expert review, although prior work has shown near-human-level performance of the same model on related neurology phenotyping tasks [24]. Future studies should include additional clinician validation of the automated phenotyping approach, ideally using independent note samples and multiple expert reviewers.

Second, all notes were drawn from a single urban academic safety-net medical center. As a result, the cohort may not be representative of MS populations seen in other clinical settings, particularly with respect to disease severity, disability burden, socioeconomic factors, access to care, and documentation practices. We also did not account for potential differences in documentation style across physicians, clinics, or health systems. The study lacked an independent internal or external validation cohort and should therefore be considered proof-of-concept rather than definitive. Future work should test whether the same four-module structure is reproducible in additional MS cohorts and across institutions.

Third, phenotype features were aggregated across all available notes for each patient to reduce variability in single-note documentation, but we did not explicitly model the number or timing of notes, nor did we examine the temporal evolution of module weights. Longitudinal analyses could determine whether module scores change over time with relapse, remission, progression, or disease-modifying therapy, and whether shifts in module composition provide clinically useful information beyond static phenotype burden.

In addition, phenotype features were scored on a binary basis as absent or present. We did not attempt to grade severity, frequency, laterality, anatomical distribution, or functional impact. As a result, a mild symptom and a severe symptom could contribute equally to the feature matrix. The feature set was intentionally limited to commonly documented MS phenotypes, and some less common manifestations, such as tinnitus, hearing loss, seizures, or other episodic symptoms, were not included. Future studies could extend this framework by incorporating severity grading, anatomical localization, temporal dynamics, and a broader set of MS-related features. As a future enhancement, phenotype features could be mapped to concepts in standard ontologies such as the Human Phenotype Ontology (HPO) or SNOMED CT.

Fourth, we focused on non-negative matrix factorization (NMF) and did not formally compare the four-module solution with alternative dimensionality-reduction, matrix-factorization, or clustering methods, such as principal components analysis, independent components analysis, factor analysis, latent class analysis, latent Dirichlet allocation/topic modeling, archetypal analysis, hierarchical or consensus clustering, Gaussian mixture models, biclustering, or graph/community-detection approaches. NMF was selected because its non-negative, parts-based representation is well suited to modeling additive symptom burden and produces clinically interpretable modules. As with all dimensionality-reduction methods, NMF necessarily involves some loss of feature-level information. In the present analysis, 17 observed phenotype features were reduced to four latent modules, and this compression was reflected in the nonzero relative reconstruction error (Table 1). Relative reconstruction error is strongly dependent on the structure of the input matrix and should not be compared directly across domains such as image matrices and sparse binary clinical phenotype matrices [36]. We therefore interpreted the four-module solution as a clinically interpretable approximation of recurring phenotype patterns rather than as a complete reconstruction of the original patient-by-feature matrix. Future studies should determine whether similar phenotype structure is recovered using alternative methods and should compare solutions using stability, reconstruction error, interpretability, and external clinical validity.

We also did not attempt to reconstruct Kurtzke Extended Disability Status Scale (EDSS) or Functional Systems Scores (FSS) from the available notes [27,28,29,37]. The notes analyzed in this study were routine clinical-care documents rather than standardized research or clinical-trial assessments, and they did not consistently contain the graded severity, functional-system scoring, ambulation distance, or assistance-level information required for valid EDSS or FSS assignment. Future studies linking note-derived module scores to prospectively collected or carefully curated EDSS and FSS measures could determine whether module dominance and phenotypic admixture provide information complementary to established MS disability scales.

Furthermore, phenotype extraction was dependent on the content and quality of clinical documentation. Physicians vary in documentation style, completeness, and attention to neurological detail, and these differences may influence the observed phenotype matrix independently of the patient’s true clinical state. Thus, some variation in phenotype burden may reflect documentation practices rather than biological or clinical differences among patients. Although physician-level documentation quality could potentially be modeled using indirect features such as note length, structure, vocabulary diversity, frequency of negative findings, or evidence of templated text, such an analysis was outside the scope of the present study.

Finally, we did not correlate module scores with independent markers of MS disease activity or progression, such as MRI lesion burden, regional atrophy, relapse rate, disability progression, serum or CSF biomarkers such as neurofilament light chain or GFAP, or genomic, immunologic, or methylomic data. Integrating NMF-derived phenotypic modules with these biological and clinical measures will be necessary to determine whether module-dominant or highly admixed phenotypes correspond to biologically meaningful MS subtypes or clinically useful prognostic groups. Future work should assess whether these modules can predict relapse, disability progression, or treatment response, or whether modules correlate with MRI, biomarker, or genetic data [25].

5. Conclusions

Multiple sclerosis lacks a single invariant phenotypic core and may be better described as a superposition of impairments across multiple functional systems. Using non-negative matrix factorization of routinely collected neurology notes, we found that MS phenotypes can be parsimoniously represented as mixtures of four latent symptom modules: sensory-visual-pain, ataxic-spastic-falls, cognitive-psychologic-fatigue, and autonomic-bladder-bowel. This representation places each patient within a two-axis interpretive framework: the dominant module identifies the patient’s predominant phenotype direction, whereas the effective number of modules quantifies the degree of phenotypic admixture across modules (Figure 3). Most patients showed admixed phenotypes, but a substantial minority had relatively pure or strongly dominant module patterns. These results provide a simple, interpretable framework for quantifying phenotypic diversity in MS and generate module-level descriptors that may be useful in downstream analyses of disease progression, MRI burden, treatment response, biomarkers, and biologically meaningful disease subtypes.