The authors propose and simulate a multistable theoretical framework. They extend the traditional low-dimensional model of smaller communities (e.g., 2 species) to a high-dimensional model of species-rich communities (e.g., 50 species). The model can account for different types of interactions, including competition and cooperation. It is impressive that the cooperation (A) and competition (B) interaction matrices can be separated and parameterized independently and randomly, which makes it convenient to test the effects of heterogeneity in different types of interactions. In particular, species are independent if both A and B matrices have non-diagonal values as zeros. It is interesting that in species-rich communities, the model with multiple stable states can replicate four ecological phenomena simply by changing the mean interaction strengths, which is a contribution to ecological theory.

After reviewing the manuscript, I recognize the importance of applying multistable dynamics to species-rich communities, a topic to some extent overlooked in ecological studies, potentially due to computational complexity. The manuscript is clearly written and accessible to readers without mathematical expertise. The theoretical analysis is coherent and systematic, and the supplementary information provides a detailed mathematical description. The illustrations are also of high quality. The identification of four distinct classes of multistability is likely to inspire additional research in this field. Nevertheless, I believe the manuscript could be further enhanced by addressing the major and several minor points below.

Major points:

1. The conclusions draw attention to the practical applications of multistability theory in predicting shifts within ecosystems. I understand that the theoretical results in this manuscript are obtained by simulations with randomly generated parameters. However, providing more information or intuition on applying the model to real-world ecosystems would be beneficial. In particular, there are several parameters that need to be quantified to implement the multistable model for a specific ecosystem. For example, what types of information are needed to parameterize the cooperative (A) and competitive (B) interaction matrices? For a pair of species, how does one disentangle the cooperative and competitive interactions that occur simultaneously between them? To what extent is it feasible to distinguish these interactions in both lab and field experiments? If it is possible, how can each type of interaction be quantified specifically?
2. If possible, I would suggest showing the applications of at least one of the four regimes as case studies.

Minor comments:

1. In Figure 2's caption, the use of 'A' and 'B' to denote both interaction matrices and figure panels may lead to confusion. It might improve readability to bold one of these notations to distinguish them clearly.
2. In lines 349-352, could you please explain in more detail how the cliques are related to ecological succession and priority effects? Is it implied that priority effects do not exist since community diversity remains consistent when transitioning among different clique states?

Thank you for your comprehensive response. I have two minor queries and suggestions for further clarification:
1. In Figure 1, I noticed a potential typographical error next to the black dot "nal stable states." Could it be that the intended phrase is "final stable states”?
2. The challenge of finding experimental data for the four distinct regimes is understandable, and I appreciate the simulations presented in Box 1 and Figure 5. Nevertheless, I was unable to find the corresponding simulation codes on GitHub. Would it be feasible for you to share these codes? As given in Box 1, it seems there is no need to parameterize matrices A and B. The only required parameter is diversity. Does this suggest that two communities with equivalent levels of diversity are anticipated to exhibit identical multistability regimes?