Background: Accuracy and precision assessment in radiomic features is important for the determination of their potential to characterize cancer lesions. In this regard, simulation of different imaging conditions using specialized phantoms is increasingly being investigated. In this study, the design and evaluation of a modular multimodality imaging phantom to simulate heterogeneous uptake and enhancement patterns for radiomics quantification in hybrid imaging is presented. Methods: A modular multimodality imaging phantom was constructed that could simulate different patterns of heterogeneous uptake and enhancement patterns in positron emission tomography (PET), single-photon emission computed tomography (SPECT), computed tomography (CT), and magnetic resonance (MR) imaging. The phantom was designed to be used as an insert in the standard NEMA-NU2 IEC body phantom casing. The entire phantom insert is composed of three segments, each containing three separately fillable compartments. The fillable compartments between segments had different sizes in order to simulate heterogeneous patterns at different spatial scales. The compartments were separately filled with different ratios of 99mTc-pertechnetate, 18F-fluorodeoxyglucose ([18F]FDG), iodine- and gadolinium-based contrast agents for SPECT, PET, CT, and T1-weighted MR imaging respectively. Image acquisition was performed using standard oncological protocols on all modalities and repeated five times for repeatability assessment. A total of 93 radiomic features were calculated. Variability was assessed by determining the coefficient of quartile variation (CQV) of the features. Comparison of feature repeatability at different modalities and spatial scales was performed using Kruskal-Wallis-, Mann-Whitney U-, one-way ANOVA- and independent t-tests. Results: Heterogeneous uptake and enhancement could be simulated on all four imaging modalities. Radiomic features in SPECT were significantly less stable than in all other modalities. Features in PET were significantly less stable than in MR and CT. A total of 20 features, particularly in the gray-level co-occurrence matrix (GLCM) and gray-level run-length matrix (GLRLM) class, were found to be relatively stable in all four modalities for all three spatial scales of heterogeneous patterns (with CQV < 10%). Conclusion: The phantom was suitable for simulating heterogeneous uptake and enhancement patterns in [18F]FDG-PET, 99mTc-SPECT, CT, and T1-weighted MR images. The results of this work indicate that the phantom might be useful for the further development and optimization of imaging protocols for radiomic quantification in hybrid imaging modalities.