To predict a protein-folding pathway, we present an alternative to the time-consuming dynamic simulation of atomistic models. We replace the actual dynamic simulation with variational optimization of a reaction path connecting known initial and final protein conformations in such a way as to maximize an estimate of the reactive flux or minimize the mean first passage time at a given temperature, referred to as MaxFlux. We solve the MaxFlux global optimization problem with an efficient graph-theoretic approach, the probabilistic roadmap method (PRM). We employed CHARMM19 and the EEF1 implicit solvation model to describe the protein solution. The effectiveness of our MaxFlux-PRM is demonstrated in our promising simulation results on the folding pathway of the engrailed homeodomain. Our MaxFlux-PRM approach provides the direct evidence to support that the previously reported intermediate state is a genuine on-pathway intermediate, and the demand of CPU power is moderate. 10.1529/biophysj.107.119214