Protein folding is a fundamental problem in modern structural biology. The nature of the problem poses challenges to the understanding of the process via computer simulations. One of the challenges in the computer simulation of proteins at the atomic level is the efficiency of sampling conformational space. Replica exchange (RE) methods are widely employed to alleviate the difficulty. To study how to best employ RE to protein folding and binding problems, We constructed a kinetic network model for RE studies of protein folding and used this simplified model to carry out "simulations of simulations" to analyze how the underlying temperature dependence of the conformational kinetics and the basic parameters of RE all interact to affect the number of folding transitions observed. When protein folding follows anti-Arrhenius kinetics, we observe a speed limit for the number of folding transitions observed at the low temperature of interest, which depends on the maximum of the harmonic mean of the folding and unfolding transition rates at high temperature. The efficiency of temperature RE was also studied on a more complicated and realistic continuous two-dimensional potential. Comparison of the efficiencies obtained using the continuous and discrete models makes it possible to identify non-Markovian effects which slow down equilibration of the RE ensemble on the more complex continuous potential. In particular, the efficiency of RE is limited by the timescale of conformational relaxation within free energy basins. The other challenges we are facing in all-atom simulations is to obtain meaningful information on the slow kinetics and pathways of folding. We present a kinetic network model which recover the kinetics using RE-generated states as the nodes of a kinetic network. Choosing the appropriate neighbors and the microscopic rates between the neighbors, the correct kinetics of the system can be recovered by running a simulation on the network.