We study the origin of robustness of yeast cell cycle cellular network through uncovering its underlying energy landscape. This is realized from the information of the steady steady probabilities by solving a discrete set of kinetic master equations for the network. We discovered that the potential landscape of yeast cell cycle network is funnelled towards the global minimum, G1 state. The ratio of the energy gap between G1 and average versus roughness of the landscape termed as robustness ratio (RR) becomes a quantitative measure of the robustness and stability for the network. The funnelled landscape is quite robust against random perturbations from the inherent wiring or connections of the network. There exists a global phase transition between the more sensitive response or less self degradation phase leading to underlying funneled global landscape with large RR, and insensitive response or more self degradation phase leading to shallower underlying landscape of the network with small RR. Furthermore, we show the more robust landscape also leads to less dissipation cost of the network. We quantify the cost in terms of the dissipation or heat loss characterized through the steady state properties: the underlying landscape and the associated flux. We found that the dissipation cost is intimately related to the stability and robustness of the network. With least dissipation cost, the network becomes most stable and robust under mutations and perturbations on sharpness of the response from input to output as well as self degradations. The least dissipation cost may provide a general design principle for the cellular network to survive from the evolution and realize the biological function. We explore the stochastic dynamics of self regulative genes from fluctuations of molecular numbers and of on and off switching of gene states due to regulatory protein binding/unbinding to the genes. We found when the binding/unbinding is relatively fast (slow) compared with the synthesis/degradation of proteins in adiabatic (non-adiabatic) case, the self regulators can exhibit one or two peak (two peak) distributions in protein concentrations. This phenomena can also be quantified through Fano factors. This shows even with the same architecture (topology of wiring), networks can have quite different functions (phenotypes), consistent with recent single molecule single gene experiments. We further found the inhibition and activation curves to be consistent with previous results in adiabatic regime, but show significantly different behaviors in non-adiabatic regimes with previous predictions with monomer binding. Such difference is due to the dimer effect and never reported before. We derive the non-equilibrium phase diagrams of mono-stability and bi-stability in adiabatic and non-adiabatic regimes. We study the dynamical trajectories of the self regulating genes on the underlying landscapes from non-adiabatic to adiabatic limit, provide a global picture of understanding and show an analogy to the electron transfer problem. We study the stability and robustness of the systems through mean first passage time (MFPT) from one peak (basin of attraction) to another and found both monotonic and non-monotonic turnover behavior from adiabatic to non-adiabatic regimes. For the first time, we explore global dissipation by entropy production and the relation with binding/unbinding processes. Our theoretical predictions for steady state peaks, fano factors, inhibition/activation curves and MFPT can be probed and tested from experiments. We established a landscape framework to explore the global stability and robustness of the dynamical systems and networks. We explore in particular a gene network motif appeared in the experimental synthetic biology studies of two genes mutually repress and activate each other with self activation and repression. We found that coherent limit cycle oscillations emerge in two regimes: adiabatic and non-adiabatic regimes, with two mechanisms of producing the stable oscillations: nonlinear cooperative interactions in the adiabatic regime and time delays due to the slow binding to the promoters in the non-adiabatic regime. In both regimes, the landscape has a topological shape of Mexican hat in protein concentrations. The shape of the Mexican hat provides the quantitative description of the capability of the system to communicate with each other in concentration space and time. Therefore, the topology of the landscape quantitatively determines the global stability and robustness of the dynamical systems and networks. The coherence of the oscillations are shown to be correlated with the shape of Mexican hat characterize by the height from the top to the ring of the hat.