Finding complex reaction and transformation pathways involving many intermediate states is, in general, not possible on the density-functional theory level with existing simulation methods, due to the very large number of required energy and force evaluations. For complex reactions, it is not possible to determine which atom in the reactant is mapped onto which atom in the product. Trying out all possible atomic index mappings is not feasible because of the factorial increase in the number of possible mappings. We use a penalty function that is invariant under index permutations to bias the potential energy surface in such a way that it obtains the characteristics of a structure seeker, whose global minimum is the reaction product. By performing a minima-hopping-based global optimization on this biased potential energy surface, we rapidly find intermediate states that lead into the global minimum and allow us to then extract entire reaction pathways. We first demonstrate for a benchmark system, namely, the Lennard-Jones cluster LJ_{38}, that our method finds intermediate states relevant to the lowest energy reaction pathway, and hence we need to consider much fewer intermediate states than previous methods to find the lowest energy reaction pathway. Finally, we apply the method to two real systems, C_{60} and C_{20}H_{20}, and show that the reaction pathways found contain valuable information on how these molecules can be synthesized.