Markov state models have become popular in the computational biochemistry and biophysics communities as a technique for identifying stationary and kinetic information of protein dynamics from molecular dynamics simulation data. In this paper, we extend the applicability of automated Markov state modeling to simulation data of molecular self-assembly and aggregation by constructing collective coordinates from molecular descriptors that are invariant to permutations of molecular indexing. Understanding molecular self-assembly is of critical importance if we want to deepen our understanding of neurodegenerative diseases where the aggregation of misfolded or disordered proteins is thought to be the main culprit. As a proof of principle, we demonstrate our Markov state model technique on simulations of the KFFE peptide, a subsequence of Alzheimer's amyloid-β peptide and one of the smallest peptides known to aggregate into amyloid fibrils in vitro. We investigate the different stages of aggregation up to tetramerization and show that the Markov state models clearly map out the different aggregation pathways. Of note is that disordered and β-sheet oligomers do not interconvert, leading to separate pathways for their formation. This suggests that amyloid aggregation of KFFE occurs via ordered aggregates from the very beginning. The code developed here is freely available as a Jupyter notebook called TICAgg, which can be used for the automated analysis of any self-assembling molecular system, protein, or otherwise.