Interstellar dust grains penetrate into the heliosphere (region in which the solar wind propagates) due to the relative motion of the Sun and the Local Interstellar Medium (LISM). Inside the heliosphere, the motion of dust particles is mainly governed by the electromagnetic force determined by the heliospheric magnetic field. Under the action of this force, the trajectories of dust grains experience intersections with each other and self-intersections. As a result, dust density accumulation regions appear. These regions are of a great interest in the context of theoretical studying and planning of upcoming space missions. The main aim of the present study is to model the interstellar dust distribution in the heliosphere and investigate the peculiarities of the number density distribution. To describe the dusty component, the cold gas model is used, while to compute the interstellar dust distribution two approaches are considered, namely, the Eulerian and Lagrangian approaches. For solving the continuity equation in the Lagrangian coordinates, the full Lagrangian method or the Osiptsov method is used. As a result, all the peculiarities of the dust distribution are investigated and it is found that they are located on caustics, i.e., the envelopes of interstellar dust trajectories. Besides it, the regular regions of overdensity (without singularities in the number density) are discovered. It is shown that the dust component accumulation regions are located in a small neighborhood of the heliospheric current sheet, at which the magnetic field changes its polarity, and in the tail of the heliosphere. The effectiveness of the Osiptsov method of solving the continuity equation is compared with the widely used Monte Carlo method (Eulerian approach). It is shown that Monte Carlo method requires extremely high resolution of computational grid to reach the level of accuracy comparable with the Osiptsov method.