Opioid use disorder is a psychological condition that affects over 200,000 people per year in the U.S., causing the Centers for Disease Control and Prevention to label the crisis as a rapidly spreading public health epidemic. It has been found that the behavioral relationship between opioid exposure and development of opioid use disorder varies greatly between individuals, implying existence of sup-populations with varying degrees of opioid vulnerability. In this study, we assessed several behavioral variables across heroin taking, refraining and seeking to establish how these factors interact with one another resulting in a heroin dependent, resilient, or vulnerable behavioral phenotype. Over 400 (male and female) heterogeneous stock rats were used in these two studies, and data were collected from two geographically distinct locations. Rats underwent heroin self-administration training, followed by a progressive ratio and heroin-primed reinstatement test. Next, rats underwent extinction training and a cue-induced reinstatement test. To assess how these variables contribute to heroin addiction vulnerability, we developed a network-based data analysis workflow. Specifically, we integrated different cohorts of rats, remove possible batch effects, and constructed a rat-rat similarity network based on their behavioral patterns. We then implemented community detection on this similarity network using a Bayesian degree-corrected stochastic block model to uncover sub-populations of rats with differing levels of opioid vulnerability. We discovered three distinct behavioral sub-populations, each with significantly different behavioral outcomes that allowed for unique characterization of each cluster in terms of vulnerability to opioid use and seeking. We implement this analysis workflow as an open source R package, named mlsbm.