The positive matrix factorization (PMF) receptor model was used for the first time to quantify the source contributions to heavy metal pollution of sediment on a national basin scale in the upstream, midstream, and downstream rivers (Teesta and Kortoya-Shitalakkah and Meghna-Rupsha and Pasur) of Bangladesh. The metal contamination status, cooccurrence, and ecotoxicological risk were also investigated. Sediment samples were collected from 30 sites at a depth range of 0 to 20 cm for analysis of 9 metals using inductively coupled plasma-mass spectrometry. The mean concentrations of metals varied for upstream, lower midstream, and downstream river segments. The results showed that chromium (Cr) exhibited a strong significant co-occurrence network with other metals (e.g., manganese [Mn], iron [Fe], and nickel [Ni]). Monte Carlo simulation results of the geo-accumulation index (Igeo; 63.3%) and risk indices (48.5%) showed that cadmium (Cd) was the main contributor to sediment pollution. However, the cumulative probabilities of sediments being polluted by metals were ranked as "moderate to heavily polluted" (Igeo 46.6%; risk index 16.7%). Toxicity unit results revealed that zinc (Zn) and Cd were the key toxic contributors to sediments. The PMF model predicted metal concentrations and identified 4 potential sources. The agricultural source (factor 1) mostly contributed to copper (Cu; 78.9%) and arsenic (As; 62.8%); Ni (96.9%) and Mn (83.5%) exhibited industrial point sources (factor 2), with 2 hot spots in northwestern and southwestern regions. Cadmium (93.5%) had anthropogenic point sources (factor 3), and Fe (64.3%) and Cr (53.5%) had a mixed source (factor 4). Spatially, similar patterns between PMF apportioning factors and predicted metal sources were identified, showing the efficiency of the model for river systems analysis. The degree of metal contamination in the river segments suggests an alarming condition for biotic components of the ecosystem.