Use the initial Center Of Mass (COM) as an origin.
It's trivial to show that the COM of a closed N-body system (in Newtonian physics) moves with a constant velocity. This is a relatively easy system to work with and transforming to other coordinate systems, such as the heliocentric one you mention, is relatively trivial compared to computing the 3D N-body system with good accuracy. A word of advice : if the math of converting between coordinate systems is a problem for you, the mathematics and algorithms of an accurate simulation will be likely much worse. I'd suggest trying out such conversion mathematics as a warm up exercise and if you're comfortable with that then dive into the simulation.
COM can be shown to move at constant velocity because it's a closed system (or can be treated as such over the timescale you mention), and so net force on the COM is zero, meaning it has to move with constant velocity.
Adding a small body like an extrasolar asteroid makes next to no difference to the system as a whole. It's a minute change. You can treat such objects as special cases, ignoring their effect on the system and just calculating the effect of the system on them. If you're planning to add very large bodies and see what happens then you need to include them in the general system all the way through the simulation.
Note there are a lot of existing simulations out there that implement a wide range of approaches. This is a problem that has been extensively studied and a search online would yield a myriad of approaches.
The motion of the Barycenter is complex. I'd avoid this complication in what is already gong to be a messy mathematical system to simulate.
Your main priorities should be gathering accurate data for initial starting values (non trivial) and choice of algorithms. Choice or ordinate system is comparatively minor.