We propose a new interpolatory framework for model reduction of large-scale bilinear systems. The input-output representation of a bilinear system in frequency domain involves a series of multivariate transfer functions, each representing a subsystem of the bilinear system. If a weighted sum of these multivariate transfer functions associated with a reduced bilinear system interpolates a weighted sum of the original multivariate transfer functions, we say that the reduced system satisfies a Volterra series interpolation. These interpolatory conditions can also ensure the necessary conditions for H2-optimal model reduction. We observe that, by carefully selecting the weights of the series, the Volterra series interpolatory conditions are transformed to the problem of interpolating a linear system with an affine parameter dependence. Such linear parametric systems can then be reduced by some method for parametric model order reduction. Linear systems where the affine parameter dependence is given as low-rank variation in the state matrix can be mapped into a non-parameterized multi-input multi-output linear system. This allows us to utilize the standard (non-parametric) linear IRKA for the problem of parameterized/bilinear interpolation. Numerical results show that the approximations are of comparable accuracy to those obtained from the bilinear iterative rational Krylov algorithm. The proposed approach, however, has the advantage that it reduces the computational costs as it involves computations associated with solving linear systems only.