A main topic in the study of topological indices is to find bounds of the indices involving several parameters and/or other indices. In this paper we perform statistical (numerical) and analytical studies of the harmonic index H(G), and other topological indices of interest, on Erd ̋os-R ́enyi (ER) graphs G(n, p) characterized by n vertices connected independently with probability p ∈ (0, 1). Particularly, in addition to H(G), we study here the (−2) sum-connectivity index χ−2(G), the modified Zagreb index M Z(G), the inverse degree index ID(G) and the Randi ́c index R(G). First, to perform the statistical study of these indices, we define the averages of the normalized indices to their maximum value: 〈H(G)〉, 〈χ−2(G)〉,〈M Z(G)〉, 〈ID(G)〉 and 〈R(G)〉. Then, from a detailed scaling analysis, we show that the averages of the normalized indices scale with the product ξ ≈ np. Moreover, we find two different behaviors. On the one hand, 〈H(G)〉 and 〈R(G)〉, as a function of the probability p, show a smooth transition from zero to n/2 as p increases from zero to one. Indeed, after scaling, it is possible to define three regimes: a regime of mostly isolated vertices when ξ 10 (H(G), R(G) ≈ n/2). On the other hand, 〈χ−2(G)〉, 〈M Z(G)〉 and 〈ID(G)〉 increase with p until approaching their maximum value, then they decrease by further increasing p. Thus, after scaling the curves corresponding to these indices display bell-like shapes in log scale, which are symmetric around ξ ≈ 1; i.e. the percolation transition point of ER graphs. Therefore, motivated by the scaling analysis, we analytically (i) obtain new relations connecting the topological indices H, χ−2, M Z, ID and R that characterize graphs which are extremal with respect to the obtained relations and (ii) apply these results in order to obtain inequalities on H, χ−2, M Z, ID and R for graphs in ER models.