The aim of this study is to develop an enhanced machine learning-based prediction models for bronchopulmonary dysplasia (BPD) and its severity through a two-stage approach integrated with the duration of respiratory support (RSd) using prenatal and early postnatal variables from a nationwide very low birth weight (VLBW) infant cohort. We included 16,384 VLBW infants admitted to the neonatal intensive care unit (NICU) of the Korean Neonatal Network (KNN), a nationwide VLBW infant registry (2013-2020). Overall, 45 prenatal and early perinatal clinical variables were selected. A multilayer perceptron (MLP)-based network analysis, which was recently introduced to predict diseases in preterm infants, was used for modeling and a stepwise approach. Additionally, we applied a complementary MLP network and established new BPD prediction models (PMbpd). The performances of the models were compared using the area under the receiver operating characteristic curve (AUROC) values. The Shapley method was used to determine the contribution of each variable. We included 11,177 VLBW infants (3,724 without BPD (BPD 0), 3,383 with mild BPD (BPD 1), 1,375 with moderate BPD (BPD 2), and 2,695 with severe BPD (BPD 3) cases). Compared to conventional machine learning (ML) models, our PMbpd and two-stage PMbpd with RSd (TS-PMbpd) model outperformed both binary (0 vs. 1,2,3; 0,1 vs. 2,3; 0,1,2 vs. 3) and each severity (0 vs. 1 vs. 2 vs. 3) prediction (AUROC = 0.895 and 0.897, 0.824 and 0.825, 0.828 and 0.823, 0.783, and 0.786, respectively). GA, birth weight, and patent ductus arteriosus (PDA) treatment were significant variables for the occurrence of BPD. Birth weight, low blood pressure, and intraventricular hemorrhage were significant for BPD ≥2, birth weight, low blood pressure, and PDA ligation for BPD ≥3. GA, birth weight, and pulmonary hypertension were the principal variables that predicted BPD severity in VLBW infants. We developed a new two-stage ML model reflecting crucial BPD indicators (RSd) and found significant clinical variables for the early prediction of BPD and its severity with high predictive accuracy. Our model can be used as an adjunctive predictive model in the practical NICU field.