  The maximum mass of neutron stars ($M_{\rm TOV}$) plays a crucial role in
understanding their equation of state (EoS). Previous studies have used the
measurements for the compactness of massive pulsars and the tidal deformability
of neutron stars in binary neutron star (BNS) mergers to constrain the EoS and
thus the $M_{\rm TOV}$. The discovery of the most massive pulsar, PSR
J0952-0607, with a mass $\sim 2.35M_{\odot}$, has provided a valuable lower
limit for $M_{\rm TOV}$. Another efficient method to constrain $M_{\rm TOV}$ is
by examining the type of central remnant formed after a BNS merger.
Gravitational wave (GW) data can provide the total mass of the system, while
accompanying electromagnetic signals can help infer the remnant type. In this
study, we combine all the previous constraints and utilize the observational
facts that about $24\%$ of the short gamma-ray bursts are followed by an X-ray
internal plateau, which indicate that roughly this fraction of BNS mergers
yield supermassive neutron stars, to perform (Markov Chain) Monte Carlo
simulations. These simulations allow us to explore the probability density
distribution of $M_{\rm TOV}$ and other parameters related to BNS mergers. Our
findings suggest that $M_{\rm TOV}$ is likely around $2.49M_{\odot} -
2.52M_{\odot}$, with an uncertainty range of approximately [$-0.16M_{\odot}$,
$0.15M_{\odot}$] ([$-0.28M_{\odot}$, $0.26M_{\odot}$]) at $1\sigma$ ($2\sigma$)
confidence level. Furthermore, we examine the type of merger remnants in
specific events like GW170817 and GW190425 to further constrain $M_{\rm TOV}$
and other relevant parameters, which can help to understand the physical
processes involved in BNS mergers.
