A numerical simulation is carried out to study the role of particles in gas-liquid-solid flows in bubble columns. An Eulerian-Lagrangian model is used and the liquid flow is modeled using a volume-averaged system of governing equations, while motions of bubbles and particles are evaluated using Lagrangian trajectory analysis. It is assumed that the bubbles remain spherical. The interactions between bubble-liquid and particle-liquid are included in the study. The discrete phase equations include drag, lift, buoyancy, and virtual mass forces. Particle-particle interactions and bubble-bubble interactions are accounted for by the hard sphere model approach. The bubble coalescence is also included in the model. Neutrally buoyant particles are used in the study. A parcel approach is used and a parcel represents a certain number of particles of same size, velocity, and other properties. Variation of particle loading is modeled by changing the corresponding number of particles in every parcel. In a previous work, the predicted results were compared with the experimental data, and good agreement was obtained. The transient flow characteristics of the three-phase flow are studied and the effects of particle loading on flow characteristics are discussed. The simulations show that the transient characteristics of the three-phase flow in a column are dominated by time-dependent vortices. The particle loading can affect the characteristics of the three-phase flows and flows with high particle loading evolve faster.