Atomistic description of protein fibril formation has remained prohibitive due to the complexity and long timescales of the conformational search problem. Here, we develop a multi-scale approach that combines a large number of atomistic molecular dynamics simulations in explicit solvent to derive Markov State Models (MSMs) for simulation of fibril growth. The search for the in-registered fully bound fibril state is modeled as a random walk on a rugged 2D energy landscape along enumerated β-sheet registry and hydrogen bonding states, whereas interconversions among nonspecific bound states and between nonspecific and hydrogen-bounded states are derived from kinetic clustering analysis. The reversible association/dissociation of an incoming peptide and overall growth kinetics are then computed from MSM trajectories. This approach is applied to derive a comprehensive description of fibril elongation of wild-type Aβ16-22 and how it is modulated by phenylalanine to cyclohexane (CHA) mutations. The resulting models recapitulate the experimental observation that mutants CHA19 and CHA1920 accelerate fibril elongation, but have a relatively minor effect on the critical concentration for fibril growth. Importantly, the kinetic consequences of mutations arise from a complex perturbation of the network of productive and non-productive pathways of fibril grown. This is consistent with the expectation that nonfunctional states will not have evolved efficient folding pathways and, therefore, will require a random search of configuration space. This study highlights the importance of describing the complete energy landscape when studying the elongation mechanism and kinetics of protein fibrils.