Neuromyelitis optica spectrum disorders lack imaging biomarkers associated with disease course and supporting prognosis. This complex and heterogeneous set of disorders affects many regions of the central nervous system, including the spinal cord and visual pathway. Here, we use graph theory-based multimodal network analysis to investigate hypothesis-free mixed networks and associations between clinical disease with neuroimaging markers in 40 aquaporin-4-immunoglobulin G antibody seropositive patients (age = 48.16 ± 14.3 years, female:male = 36:4) and 31 healthy controls (age = 45.92 ± 13.3 years, female:male = 24:7). Magnetic resonance imaging measures included total brain and deep grey matter volumes, cortical thickness and spinal cord atrophy. Optical coherence tomography measures of the retina and clinical measures comprised of clinical attack types and expanded disability status scale were also utilized. For multimodal network analysis, all measures were introduced as nodes and tested for directed connectivity from clinical attack types and disease duration to systematic imaging and clinical disability measures. Analysis of variance, with group interactions, gave weights and significance for each nodal association (hyperedges). Connectivity matrices from 80% and 95% F-distribution networks were analyzed and revealed the number of combined attack types and disease duration as the most connected nodes, directly affecting changes in several regions of the central nervous system. Subsequent multivariable regression models, including interaction effects with clinical parameters, identified associations between decreased nucleus accumbens (β = −0.85, P = 0.021) and caudate nucleus (β = −0.61, P = 0.011) volumes with higher combined attack type count and longer disease duration, respectively. We also confirmed previously reported associations between spinal cord atrophy with increased number of clinical myelitis attacks. Age was the most important factor associated with normalized brain volume, pallidum volume, cortical thickness and the expanded disability status scale score. The identified imaging biomarker candidates warrant further investigation in larger-scale studies. Graph theory-based multimodal networks allow for connectivity and interaction analysis, where this method may be applied in other complex heterogeneous disease investigations with different outcome measures.