This study investigates an integrated train operation plan problem with flexible train compositions and multiple service routes, which allows trains to couple, decouple and turn around at intermediate turnaround stations based on virtual coupling technology. An integrated optimization model for demand-oriented train timetables and rolling stock circulation plans is developed with the consideration of train service routes and the number of train services, by minimizing both passenger waiting time and operation costs. Specifically, with the path-based space-time network representation, we formulate a novel rolling stock circulation model, which accounts for coupling and decoupling operations and stabling track capacity constraints with lower complexity than the assignment model. In addition, we improve the linear dynamic passenger flow loading model, which extends the operational application with multiple service routes and characterizes various passenger waiting behaviors. To solve the proposed model, an adaptive simulated annealing (ASA) algorithm is designed to obtain high-quality solutions using flow-oriented and random operators. Lastly, the performance of the proposed models and algorithms is verified by small-scale numerical experiments. Then, the efficiency of the proposed approaches is further demonstrated through a real-world case study based on Beijing Subway Changping Line, showing a 42% reduction in total passenger waiting time, a decrease of 16 rolling stock, and a 10% reduction in variable operation costs compared to the current operation mode.