Recently, the fractional Fokker-Planck equations (FFPEs) with multiple internal states are built for the particles undergoing anomalous diffusion with different waiting time distributions for different internal states, which describe the distribution of positions of the particles [Xu and Deng, Math. Model. Nat. Phenom., $\mathbf{13}$, 10 (2018)]. In this paper, we first develop the Sobolev regularity of the FFPEs with two internal states, including the homogeneous problem with smooth and nonsmooth initial values and the inhomogeneous problem with vanishing initial value, and then we design the numerical scheme for the system of fractional partial differential equations based on the finite element method for the space derivatives and convolution quadrature for the time fractional derivatives. The optimal error estimates of the scheme under the above three different conditions are provided for both space semidiscrete and fully discrete schemes. Finally, one- and two-dimensional numerical experiments are performed to confirm our theoretical analysis and the predicted convergence order.