The accurate representation of the ionospheric electron density in 3-dimensions is a challenging problem because of the nature of horizontal and vertical structures on both small and large scales. This paper presents the development of a global three-dimensional (3-D) electron density reconstruction based on radio occultation data during 2006–2019 and neural networks. We demonstrate that the developed model based on COSMIC dataset only is capable of reproducing different ionospheric features when compared to independent datasets from ionosondes and incoherent scatter radars (ISR) in low, middle and high latitude regions. Following some existing modelling efforts based on similar or related datasets and technique we divided the problem into fine resolution grid cells of 5∘×15∘ (geographic latitudes/longitudes) followed by development of the neural network subroutine per cell and later combining all the 864 sub-models to compile one global model. This approach has been demonstrated to be appropriate in enabling neural networks to learn, reproduce and generalise local and global behaviour of the ionospheric electron density. Based on ISR data, the 3D model improves maximum electron density of the F2 layer (NmF2) prediction by 10%–20% compared to IRI 2016 model during quiet conditions. For estimation of ionosonde ordinary critical frequency of the F2 layer (foF2) in 2009 at 1200 UT (universal time), the developed 3-D model gives average root mean square error (RMSE) values of 0.83 MHz, 1.06 MHz and 1.16 MHz compared to the IRI 2016 values of 0.92 MHz, 1.09 MHz and 1.01 MHz over the Africa–European, American and Asian sectors respectively making their performances statistically comparable. Compared to ionosonde data, the IRI 2016 model consistently shows a better performance for the hmF2 modelling results in almost all sectors during the investigated periods.