Supernovae (SNe), the most energetic stellar feedback mechanism, are crucial for regulating the interstellar medium (ISM) and launching galactic winds. We explore how supernova remnants (SNRs) create a multiphase medium by performing three-dimentional hydrodynamical simulations at various SN rates, S, and ISM average densities, (n) over bar. The evolution of an SNR in a self-consistently generated three-phase ISM is qualitatively different from that in a uniform or a two-phase warm/cold medium. By traveling faster and further in the low-density hot phase, the domain of an SNR increases by > 10(2.5). Varying (n) over bar and S, we find that a steady state can only be achieved when the hot gas volume fraction f(V,hot) less than or similar to 0.6 +/- 0.1. Above that level, overlapping SNRs render connecting topology of the hot gas, and the ISM is subjected to thermal runaway. Photoelectric heating (PEH) has a surprisingly strong impact on f(V,hot). For (n) over bar greater than or similar to 3 cm(-3), a reasonable PEH rate is able to suppress the thermal runaway. Overall, we determine the critical SN rate for the onset of thermal runaway to be S-crit = 200((n) over bar /1 cm(-3))(k) (E-SN/10(51) erg)(-1) kpc(-3) Myr(-1), where k = (1.2, 2.7) for (n) over bar <= 1 and > 1 cm(-3), respectively. We present a fitting formula of the ISM pressure P((n) over bar ,S), which can be used as an effective equation of state in cosmological simulations. Despite the five orders of magnitude span of ((n) over bar, S), the average Mach number varies little: M approximate to 0.5 +/- 0.2, 1.2 +/- 0.3, and 2.3 +/- 0.9 for the hot, warm, and cold phases, respectively.