Both atomistic and mesoscale simulation techniques have been extensively employed to gain fundamental understanding of the structures, deformation mechanisms, and structure-property relationships in bulk metallic glasses (BMGs), each with its unique strengths and limitations. Nevertheless, there is a limited degree of synergistic integration between the two approaches. In this study, we extract key properties of shear transformation zones (STZs) directly from the atomistic simulations, including their size, number of shear modes, eigenstrain, and most importantly, the activation energy barrier spectrum as a function of cooling history and strain rate. We then incorporate these STZ properties into a heterogeneously randomized STZ dynamic model implemented in a kinetic Monte Carlo algorithm to study parametrically the deformation microstructure, shear band formation and stress-strain behavior of BMGs. Two important characteristics of STZ activation that dictate the strength and ductility of a glass are identified. One is the average of the activation energy barrier spectrum (approximated by a Gaussian distribution), determined by the glass composition and processing history such as the cooling rate. The other is the amount of shift of the Gaussian distribution towards smaller activation energy barrier values during deformation, which is determined by the initial structural states and strain rate during deformation, and exhibits a saturation value. These findings have allowed us to gain important fundamental insights into the correlation between the degree of shear-induced softening and the general deformation behavior of BMGs, leading to a better understanding of the correlation between the processing history/loading condition and the mechanical behavior.