We will discuss a parallel boundary integral method for simulating highly concentrated vesicle suspensions in a Stokesian fluid. The simulation of high volume fraction vesicle suspensions which are representative of real biological systems (such as blood with 35% - 50% volume fraction for RBCs) presents several challenges. It requires computing accurate vesicle-vesicle interactions at length scales where standard quadratures are too expensive. The inter-vesicle separation can become arbitrarily small leading to vesicle collisions. Numerical errors can accumulate over time, making long-timescale simulations inaccurate. We tackle these challenges by developing state-of-the-art parallel algorithms for efficient computation of boundary integrals and an adaptive time-stepping scheme. We will also discuss algorithms for handling vesicle collisions and re-meshing of vesicle surfaces.