Monolithic bulk metallic glasses (BMGs) exhibit a unique combination of mechanical properties, such as high strength and large elasticity limits, but the lack of ductility is considered the main Achilles' heel of BMG systems. To increase the competitiveness of BMGs vis-à-vis conventional structural materials, the problem of catastrophic failure via intense plastic strain localization (‘shear banding’) has to be addressed. Recent experimental observations suggest that the addition of structural heterogeneities, in the form of crystalline particles, to BMG systems hinders the catastrophic propagation of shear bands and leads to enhanced ductility. These structural heterogeneities can be introduced by either forming BMG composites, where second-phase crystalline particles accommodate applied loads via martensitic transformation mechanisms, or developing glassy alloys that precipitate crystalline particles under deformation, a process by which further deformation can be sustained by twinning mechanisms in the crystalline phase. In this work, we present a non-linear continuum model capable of capturing the structural heterogeneity in the glassy phase and accounting for intrinsic work hardening via martensitic transformations in second-phase reinforcements in BMG composites and deformation twinning in precipitated crystalline particles. Simulation results reveal that in addition to intrinsic work hardening in the crystalline phase, particle size greatly affects the overall mechanical behavior of these BMG systems. The precipitation of crystalline particles in monolithic BMGs yields two-phase microstructures that promote more homogeneous deformation, delay the propagation of incipient shear bands, and ultimately result in improved ductility characteristics.