We present a class of exponential integrators to compute solutions of the stochastic Schrödinger equations arising from the modeling of open quantum systems. To be able to implement the methods within the same framework as the deterministic counterpart, we express the solution using Kunita's representation. With appropriate truncations, the solution operator can be written as matrix exponentials, which can be efficiently implemented by the Krylov subspace projection. The accuracy is examined in terms of the strong convergence by comparing trajectories, and in terms of the weak convergence by comparing the density-matrix operators. We show that the local accuracy can be further improved by introducing third-order commutators in the exponential. The effectiveness of the proposed methods is tested using the example from Di Ventra et al. [J. Phys.: Condens. Matter 16, 8025 (2004)JCOMEL0953-898410.1088/0953-8984/16/45/024].