We study the inertial migration of finite-size neutrally buoyant spherical particles in dilute and semi-dilute suspensions in laminar square duct flow. We perform several direct numerical simulations using an immersed boundary method to investigate the effects of the bulk Reynolds number Re b , particle Reynolds number Re p and duct to particle size ratio h/a at different solid volume fractions φ, from very dilute conditions to 20%. We show that the bulk Reynolds number Re b is the key parameter in inertial migration of particles in dilute suspensions. At low solid volume fraction (φ = 0.4%) and low bulk Reynolds number (Re b = 144), particles accumulate at the center of the duct walls. As Re b is increased, the focusing position moves progressively towards the corners of the duct. At higher volume fractions, φ = 5, 10 and 20%, and in wider ducts with Re b = 550, particles are found to migrate away from the duct core towards the walls. In particular, for φ = 5 and 10%, particles accumulate preferentially at the corners. At the highest volume fraction considered, φ = 20%, particles sample all the volume of the duct, with a lower concentration at the duct core.For all cases, we find that particles reside longer times at the corners than at the wall centers. In a duct with lower duct to particle size ratio h/a (i.e. with larger particles), Re b = 144 and φ = 5% we find that particles preferentially accumulate around the corners. Hence, the volume fraction plays a key role in defining the final distribution of particles in semi-dilute suspensions. The presence of particles induces secondary cross-stream motions in the duct cross-section, for all φ. The intensity of these secondary flows depends strongly on particle rotation rate, on the maximum concentration of particles in focusing positions, and on the solid volume fraction. We find that the secondary flow intensity increases with the volume fraction up to φ = 5%. However, beyond φ = 5% excluded volume effects lead to a strong reduction of cross-stream velocities. Inhibiting particles from rotating also results in a substantial reduction of the secondary flow intensity, and in variations of the exact location of the focusing positions.