Abstract: In voxel-based neuroimaging disease prediction, it was recently found that in addition to lesion features, there exists another type of feature called "Procedural Bias", which is introduced during preprocessing and can further improve the prediction power. However, traditional sparse learning methods fail to simultaneously capture both types of features due to their heterogeneity in sparsity types. Specifically, the lesion features are spatially coherent and suffer from volumetric degeneration, while the procedural bias refers to enlarged voxels that are dispersedly distributed. In this paper, we propose a new method based on differential inclusion, which generates a sparse regularized solution path on a couple of parameters that are enforced with heterogeneous sparsity to capture lesion features and the procedural bias separately. Specifically, we employ Total Variation with a non-negative constraint for the parameter associated with degenerated and spatially coherent lesions; on the other hand, we impose $\ell_1$ sparsity with a non-positive constraint on the parameter related to enlarged and scatterly distributed procedural bias. We theoretically show that our method enjoys model selection consistency and $\ell_2$ consistency in estimation. The utility of our method is demonstrated by improved prediction power and interpretability in the early prediction of Alzheimer's Disease.