The work is devoted to the development of a fundamentally new way of modeling the ionospheric D-region – deterministic-probabilistic. The results of Ne calculations using this technique are analyzed. Research of this kind is of fundamental importance, related to the rejection of a purely deterministic description of a continuously changing environment such as the ionosphere. In this work, the electron density is calculated using a five-component system of ionization-recombination cycle equations. Probability density functions (PDFs) of input parameters of the model are used to solve the system. The most important sources of the D-region ionization are taken into account to calculate PDFs of the ionization rate. The necessary number of iterations is determined by the convergence of PDFs of the electron density from 50 km to 85 km at midlatitudes under different heliogeophysical conditions. Theoretical Ne PDFs have been shown to be in good agreement with two experimental databases on electron density, especially at large D-region heights. The next important stage of modeling is the thorough verification of Ne PDFs from experimental radiophysical data on VLF–LF propagation.