Metallic glasses (MGs) possess remarkably high strength but often display only minimal tensile ductility due to the formation of catastrophic shear bands. Purposely enhancing the inherent heterogeneity to promote distributed flow offers new possibilities in improving the ductility of monolithic MGs. Here, we report the effect of the spatial heterogeneity of elasticity, resulting from the inherently inhomogeneous amorphous structures, on the deformation behavior of MGs, specifically focusing on the ductility using multiscale modeling methods. A highly heterogeneous, Gaussian-type shear modulus distribution at the nanoscale is revealed by atomistic simulations in Cu64Zr36 MGs, in which the soft population of the distribution exhibits a marked propensity to undergo the inelastic shear transformation. By employing a mesoscale shear transformation zone dynamics model, we find that the organization of such nanometer-scale shear transformation events into shear-band patterns is dependent on the spatial heterogeneity of the local shear moduli. A critical spatial correlation length of elastic heterogeneity is identified for the simulated MGs to achieve the best tensile ductility, which is associated with a transition of shear-band formation mechanisms, from stress-dictated nucleation and growth to structure-dictated strain percolation, as well as a saturation of elastically soft sites participating in the plastic flow. This discovery is important for the fundamental understanding of the role of spatial heterogeneity in influencing the deformation behavior of MGs. We believe that this can facilitate the design and development of new ductile monolithic MGs by a process of tuning the inherent heterogeneity to achieve enhanced ductility in these high-strength metallic alloys.